# Mode Tracking¶

The Mode Tracking task is an interface for mode- and intensity-tracking methods, adapted from the MoViPac suite [6- 10]. These methods can be used to obtain select normal modes, without having to calculate the entire vibrational spectrum. It does this through an iterative procedure.

Warning

There was a bug in the Mode Tracking implementation in early versions of AMS2019, leading to unreliable convergence. The problem has been fixed in the AMS2019.103 subrelease. We advise users not to run Mode Tracking calculations with the AMS2019.101 and AMS2019.102 subreleases.

Mode Tracking starts with a numerical frequency calculation, which refines the initial guess $$b^m$$ for the selected mode. The error of this mode with respect to the true Hessian eigenvector is calculated. This error is used in a (Jacobi-)Davidson algorithm to generate an additional mode. In subsequent iterations, we use these modes as approximations to the true normal modes. In this way, the error of the mode is minimised iteratively, yielding a closer approximation to true normal modes. This is how Mode Tracking differs from the Mode Refinement methods, in that it guarantees that the obtained modes are correct (assuming the procedure has converged).

The GUI tutorial on Mode Tracking.

## Theory¶

During the numerical frequency calculation, we obtain also the projection of the Hessian onto this mode:

$\sigma_i = H^m \cdot b_i^m = \frac{\partial{}^{2}E}{\partial{}R_i^{m}\partial{}b^m}$

This term is calculated through finite differences on the analytical gradients of the electronic energy along the mass-weighted normal modes $$q^m$$. The index $$i$$ denotes the $$3N$$ nuclear coordinates. From this projection a Rayleigh matrix is generated:

$\tilde{H}^m = {B^m}^T \cdot \Sigma$

Here, $$B^m$$ and $$\Sigma$$ are matrices containing the $$b^m$$ and $$\sigma$$ vectors for all foregoing iterations. During each iteration $$k$$, if we have not converged, we generate an updated guess vector $$b_k^m$$, and so the number of vectors in the matrices above is equal to the number of iterations $$k$$. The eigenvectors of $$\tilde{H}^m$$ give us the coefficient series for linear combinations of the guess modes $$b^m$$ such that we obtain approximations for the true normal modes:

$Q^m = \sum_k c_k \cdot b_k^m$

Each iteration, we expand the vector basis $$B^m$$, which allows this series expansion to come closer to the true normal modes each time. We can also calculate the error of this mode with respect to how close it is to being an eigenvalue of the real Hessian:

$r = \sum_k c_k \cdot \Big[\sigma_k - \lambda \cdot b_k\Big]$

Here, $$\lambda$$ is the corresponding eigenvalue of $$\tilde{H}^m$$. $$r$$ is the residual vector, giving the error for each vector element. It should be zero if the mode is an exact eigenvector of the true Hessian.

Since $$\tilde{H}^m$$ may give multiple eigenvectors, several approximate modes will be generated during those iterations. Out of these, 1 mode is identified as the mode of interest according to the specified tracking method. If the residual of this mode has been minimised sufficiently, the procedure has converged. If not, we generate a new guess vector $$b_k^m$$. There are 2 algorithms for generating this new guess, set by UpdateMethod:

### Davidson method:¶

The Davidson method uses a pre-conditioner $$D$$ to generate a new guess mode from the residual vector of the mode selected by the tracking method:

$b_k^m = D^{-1} \cdot r$

This preconditioner is constructed from an approximation of the Hessian:

$D = H_A - \lambda \cdot I$

The Davidson method works reasonably well, but can have trouble converging if the approximate modes or the Hessian are too accurate. This results as the new vectors that are generated do not necessarily extend the span of the basis. 

### vdVorst-Sleijpen-Jacobi-Davidson:¶

This variant of the Jacobi-Davidson scheme from Sleijpen & vdVorst  automatically makes the new guess vector orthogonal to the normal mode selected by the tracking method:

$b_k^m = \Big(\frac{Q^{m}D^{-1}r}{Q^{m}D^{-1}Q^m}\Big)D^{-1}Q^m - D^{-1}r$

The new vector is therefore guaranteed to extend the span of the basis as much as possible, and thus also eliminates the aforementioned issue with the Davidson method. In general, it is therefore recommended to use this Jacobi-Davidson method since it is found to converge faster, and be more reliable, as a result of yielding better guess modes.

There are 4 methods to obtain the approximate Hessian $$H_A$$, used by both update methods. They are set by HessianGuess:

• UFF is the default, which generates the approximate Hessian using UFF. While this Hessian may not yield the correct modes by itself, it produces good results as a preconditioner since it correctly represents the molecular structure.

• File will read the Hessian from an AMS output file, which can be specified in HessianPath. Using a Hessian from a more advanced method will generally yield better results for the Jacobi-Davidson method. The Davidson method will however experience difficulties with convergence as the Hessian becomes too accurate. 

• Inline will read a Hessian specified in the input file, in the HessianInline block. This allows you to use Hessians generated in external programs:

ModeTracking

HessianGuess Inline

# Approximate Hessian for H2O: 3 x nAtoms = 9 so 9x9 Hessian

HessianInline
0.62088786   0.00000000   0.00000000  -0.31044393   0.00000000  -0.21902068  -0.31044393   0.00000000   0.21902068
0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000
0.00000000   0.00000000   0.32143213  -0.15284008   0.00000000  -0.16071607   0.15284008   0.00000000  -0.16071607
-0.31044393   0.00000000  -0.15284008   0.33598889   0.00000000   0.18593038  -0.02554496   0.00000000  -0.03309030
0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000
-0.21902068   0.00000000  -0.16071607   0.18593038   0.00000000   0.15761846   0.03309030   0.00000000   0.00309761
-0.31044393   0.00000000   0.15284008  -0.02554496   0.00000000   0.03309030   0.33598889   0.00000000  -0.18593038
0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000
0.21902068   0.00000000  -0.16071607  -0.03309030   0.00000000   0.00309761  -0.18593038   0.00000000   0.15761846
End

End

• Unit uses the unit matrix. This is evidently not a good approximation for the Hessian, and is not intended to be used for proper Mode Tracking runs. However: using a poor approximation for the Hessian can result in basis vectors being generated that we would not obtain otherwise. Running Mode Tracking with this option can allow you to “probe” the vector space to obtain guesses for normal modes, which can be used as starting points for proper Mode Tracking calculations. It is however generally recommended to instead do e.g. a DFTB or UFF run if your goal is to obtain guess modes.

In later iterations, the basis $$B^m$$ will become larger. In order to improve the guess modes even further, an iterative Gram-Schmidt procedure is used to orthogonalise the new guess mode to the existing basis. An iterative procedure is necessary to account for numerical noise.

• GramSchmidt sets whether to perform this Gram-Schmidt orthogonalisation step. It is on by default.
• GramSchmidtTolerance is the absolute tolerance for orthogonality of the basis. It is evaluated with respect to the norm of the overlap vector between the new guess mode and the basis of the previous iteration $$||{b_k^m}^T B^m||$$.
• GramSchmidtIterations is the maximum number of allowed iterations during the Gram-Schmidt procedure.

The default settings for the Gram-Schmidt procedure should work for almost all systems.

• Displacement is the displacement stepsize (in Bohr) that is used for calculating frequencies, IR intensities and the Hessian projections through finite differences. The stepsize is provided for displacements along the Cartesian normal modes.
• MaxIterations is the maximum allowed number of iterations that the Mode Tracking calculation may go through. If this number is reached, the calculation will stop even if convergence was not achieved. If no value is supplied, a default of $$3N/2$$ will be used. This is approximately the maximum number of iterations where the procedure remains computationally competitive with the construction of the full Hessian.
• The ModeScanning key can be used to automatically run a numerical frequencies calculation on the new modes $$Q$$ after the Mode Tracking calculation has finished. Ritz vectors are obtained here as linear combinations of the guess modes, which in turn follow from finite difference calculations. This makes it possible for numerical errors to accumulate in the normal modes. Only in exceptional cases will these errors be significant, and running this additional refinement step is therefore only necessary if you need complete certainty that the results are accurate.

## Tracking methods¶

The TrackingMethod parameter allows you to select what property of the normal modes you want to track. At the end of each iteration, we obtain a set of approximate normal modes. The tracking method identifies which of these modes fits best for some criterium, and either returns this mode as the calculation result, or, if convergence was not achieved, uses it to generating a new basis mode for the next iteration. In general these methods are distinguished in 2 categories:

### Mode Tracking:¶

The original tracking methods focus on obtaining as accurate as possible a normal mode for the system. This class of tracking methods focusses either on accuracy of the mode, or obtaining modes with particular vibrational character:

• OverlapInitial is the default tracking method. Here, we choose the mode which resembles most closely the guess mode that was initially supplied $$b_1^m$$. This is done by choosing the mode which has the greatest overlap with the initial guess vector. This method allows us to direct the optimisation towards modes that e.g. involve particular atoms or include particular bending/stretching vibrations.
• OverlapPrevious instead chooses the mode which resembles closest the approximate mode of the previous iteration $$Q_k^m$$. This procedure allows a bit more flexibility in the optimisation. Since we essentially “forget” about earlier iterations, this procedure allows the optimisation to correct errors in the initial guess. (It is possible for instance that the initial guess included 2 different bond stretches which do not mutually occur in the true modes. This method will then converge quicker to a mode involving only 1 of these stretches, whereas OverlapInitial will take a much larger number of iterations to achieve this, if it does so at all.) Do note that this means that the final mode that you obtain does not necessarily represent the mode you initially supplied.
• DifferenceInitial works the same as OverlapInitial, except that it chooses the mode which has the smallest norm for the difference vector between the initial mode and the approximate normal modes of this iteration. The use of the difference vector prioritises deviations in the dominant parts of the vibrational character. E.g. if a mode consists primarily of a CO stretch, plus some minor vibrations in a carbon backbone, it may be desired to prioritise getting the correct force constant for the dominant CO stretch. This is achieved using these difference vector methods. In general, overlap methods still work well in these situations, and the use of difference methods should only be necessary in extreme cases.
• DifferencePrevious is also the same as DifferenceInitial except for the use of the difference vector norm as the selection criterium.
• FreqInitial chooses the mode with the frequency closest to that of the initial guess. This allows us to direct the tracking towards modes in a particular frequency region of the spectrum. Note that convergence for these frequency-based methods is slightly slower since the character of the mode itself is not included in the selection criteria, allowing for larger differences in the modes between iterations.
• FreqPrevious is similar to FreqInitial except that we choose the mode with the frequency closest to that of the previous iteration. This allows the optimisation more freedom to move away from the frequency region of the initial guess, and thus allows to correct somewhat for poor initial guesses.
• HighestFreq chooses the mode with the highest frequency. This method can be used if it is desired to track particular characteristic high-frequency vibrations.
• LowestResidual chooses the mode which has the smallest norm for the residual vector (see the ‘Convergence’ section below.) This method only focusses on obtain the most accurate mode, regardless of vibrational character or where it lies in the spectrum. This method should generally only be used as a pre-conditioner if you have very little information on what the normal modes should look like. (Since it is basically a non-directed optimisation.) This method will then try and find the normal mode closest to your guess. The approximate normal mode obtained this way will most likely not have converged yet, but should give you an indication of what the normal modes may look like. You can use these modes to refine your initial guess, and then do a new Mode Tracking run using any of the other tracking parameters to obtain the desired mode. Although this strategy is possible, it is generally recommended to use an approximate method to get an initial guess for the normal modes instead (as shown in the examples).

### Intensity Tracking:¶

This class of methods focusses on tracking modes based on their intensity in e.g. the infrared spectrum, rather than focussing on getting a mode with a particular type of vibration.

• IRInitial chooses the mode with the IR intensity closest to that of the initial guess. This constrains the optimisation to modes which are IR active, a property that may be lost when using mode tracking update methods.
• IRPrevious similarly chooses the mode with the IR intensity closest to that of the previous iteration. This allows the method some more flexibility in varying the intensity of the vibration, and thus works better if the initial guess is not that good.
• HighestIR chooses the mode with the highest IR intensity. This option can be used to find the modes associated with sharp peaks in the IR spectrum.

With Intensity Tracking, we essentially add an additional requirement to the modes: they must have a particular IR intensity. This constrained search has different convergence characteristics than conventional mode tracking, which you should take into account when setting up the mode tracking calculations.

• The majority of modes will have near-zero IR intensity. If we use a near-zero IR intensity mode as our initial guess, and request IRIntitial or IRPrevious, then we could be tracking any of one of these. Conversely, convergence behaviour will be poor since the generated basis modes are essentially random. If you are trying to obtain a high IR-intensity mode, use an IR-susceptible mode.

Note

In our conventional work-flow, we recommend starting mode tracking or refinement calculations from a set of approximate normal modes obtained from a semi-empirical or force-field method. Note however, that these method often do not produce accurate IR intensities. When selecting the initial guess mode, do not use the IRRange or related options in the ModeSelect block. This will cause you to miss vibrations which were incorrectly labelled with low IR intensity, or vice versa. Instead, rely on chemical intuition to identify the modes which contain commonly IR active vibrational components (such as C-O or N-H stretches). You can use ADFSpectra in the GUI to visualise the vibrational modes, to help you in this process.

• To allow the intensity tracking procedure to converge faster, it is recommended to use the IRPrevious tag instead of the IRInitial tag. As discussed earlier, the former allows more flexibility in the optimisation procedure, which counters the rigidity imposed by the intensity constraint. Intensity tracking methods often need this additional flexibility in generating guess modes to converge to the desired modes.
• Poor Initial Guesses: During each iteration, we still use the mode tracking methods to generate new basis modes. These basis modes try to expand the span of the basis with respect to the vibrational character of the modes. Note that this expansion does not guarantee that we will expand the basis specifically in the sub-span of IR-susceptible vibrations. If the initial guess for intensity tracking is correct, we already start our search in the sub-span vicinity of the normal modes. Basis expansion is then more efficient and there is a high chance that new guess modes sample the IR characteristic vibrations. For intensity tracking it is therefore discouraged to use poor initial guess modes.
• HighestIR is considered a “pure” intensity tracking method, in that it is used specifically to target characteristics of the IR spectrum irrespective of the underlying vibrational character. Consequently, the normal mode character can vary a lot between iterations. In order to assure that the procedure converges to the desired modes, it is recommended to use sufficiently strict tolerances (see the Convergence section). If the tolerances are too lax, the program may consider the modes to be “good enough” based on residual minimisation, even though there may be another mode with a higher IR intensity. For this reason it is generally recommended to use ToleranceForNorm values 1 order of magnitude lower than the default, or around 0.00005.

## Selecting modes¶

It is possible to track multiple modes in a single Mode Tracking calculation. The Mode Tracking task will then run the Mode Tracking algorithm for each mode in order.

The initial guess for the mode which will be tracked can be supplied in several ways. This is governed by TrackedMode:

• Inline will make the module read the mode from the input file. If this option is selected, you can supply the mode in the ModeInline block. It is possible to supply multiple modes by adding additional ModeInline blocks. The modes are given with one line for the x,y,z-displacement per atom, and in the same order, as the Atoms block in System:

ModeTracking

TrackedMode Inline

ModeInline
0.00000000   0.00000000   -0.03815965
-0.18888544   0.00000000    0.30281066
0.18888544   0.00000000    0.30281066
End

ModeInline
0.00000000   0.00000000   -0.02243153
0.32132452   0.00000000    0.17800237
-0.32132452   0.00000000    0.17800237
End

...

End

• File will make the module read modes from an AMS or engine output file, specified by ModePath. Modes generated using DFTB can be read from the dftb.rkf file and optimised using Mode Tracking for example. When this option is selected, all the vibrational modes present in the file are read first. The ModeSelect block then specifies for which of these modes you want to perform the Mode Tracking calculation.

• Hessian will generate modes as the eigenvectors of the approximate Hessian selected for the preconditioner in HessianGuess. This also allows modes to be generated for Hessians obtained from external programs. ModeSelect specifies which of the generated vibrational modes are selected for Mode Tracking.

• Settings for the ModeSelect block are discussed on the main page.

• MassWeighInlineMode decides whether the initial guess modes need to be mass-weighed. As discussed above, Mode Tracking uses mass-weighted normal modes. In most cases, the normal modes are given in regular Cartesian coordinates however. By setting MassWeighInlineMode true, these Cartesian modes are converted into mass-weighted modes by the program. If you supply a mass-weighted mode through the ModeInline block however, you do not need the program to do the mass-weighing, and you should set MassWeighInlineMode false.

## Convergence¶

In order to guide the Mode Tracking procedure, several convergence criteria are used:

• ToleranceForNorm is the absolute tolerance for convergence of the norm of the residual vector. The residual vector is a vector containing the error for each element of the normal mode, and we use the norm as a measure for the total error. If the total error is smaller than this threshold, we consider the mode to be a true normal mode and we stop iterating. Since the value of this norm depends on the length of the residual vector hence the number of atoms in the system, this tolerance is scaled internally to the number of atoms. 0.0005 is used as a default value for which most systems will converge to reasonably accurate modes in not too many iterations. If you want a more accurate approximation, you can decrease this value by e.g. 1 order of magnitude. (Consider running using the default settings, and reading the norm at convergence from the logfile. The new norm can be chosen to be lower than this value to ‘force’ the method into another iteration.)
• ToleranceForResidual is the absolute tolerance for the maximum component of the residual vector. Particularly in larger systems, where the vibration may be dominated by a small number of atoms, the error associated with the vibration of the majority of atoms may be small (the scaled residual norm will be small). The error for the atoms involved in the vibration may be comparatively large then, which is why we also check convergence for the maximum component of the error. Note that both the norm and this max. error are checked simultaneously. By varying strictness of the criteria for the norm and the max. error separately, you can prioritise either the total vibration or more localised character.
• ToleranceForBasis checks that the basis mode generated in the previous iteration, through the (Jacobi-)Davidson method, contributes to the approximate normal mode. Since the approximate mode is taken as a linear combination of the basis modes, its linear expansion coefficient must be larger than this tolerance.

Note that the iterative procedure is stopped as soon as any one of these convergence criteria is satisfied. The default values for these parameters should be applicable for most cases, but can be adjusted if needed. If stricter criteria are required, it is recommended to adjust both ToleranceForNorm and ToleranceForResidual.

## Overview of input options¶

Below is the overview for all the keys in the ModeTracking block:

ModeTracking
Displacement float
HessianGuess [Unit | File | UFF | Inline]
HessianInline # Non-standard block. See details.
...
End
HessianPath string
MassWeighInlineMode [True | False]
MaxIterations integer
ModeInline # Non-standard block. See details.
...
End
ModePath string
ModeSelect
FreqAndIRRange float_list
FreqRange float_list
Full [True | False]
HighFreq integer
HighIR integer
IRRange float_list
ImFreq [True | False]
LowFreq integer
LowFreqNoIm integer
LowIR integer
ModeNumber integer_list
End
ScanModes [True | False]
ToleranceForBasis float
ToleranceForNorm float
ToleranceForResidual float
TrackedMode [Inline | File | Hessian]
TrackingMethod [...]
UpdateMethod [JD | D]
End

ModeTracking
Type: Block Input data for ModeTracking task.
Displacement
Type: Float 0.01 Step size (in Bohr) for finite difference calculation of frequencies and IR intensities during mode tracking iterations.
HessianGuess
Type: Multiple Choice UFF [Unit, File, UFF, Inline] Sets how to obtain the guess for the Hessian used in the preconditioner.
HessianInline
Type: Non-standard block Initial guess for the (non-mass-weighted) Hessian in a 3N x 3N block, used when [HessianGuess] = [Inline].
HessianPath
Type: String Path to a .rkf file containing the initial guess for the Hessian, used when [HessianGuess] = [File].
MassWeighInlineMode
Type: Bool True The supplied modes must be mass-weighed. This tells the program to mass-weigh the supplied modes in case this has not yet been done. (True means the supplied modes will be mass-weighed by the program, e.g. the supplied modes are non-mass-weighed.)
MaxIterations
Type: Integer Maximum number of allowed iterations.
ModeInline
Type: Non-standard block True Coordinates of the mode which will be tracked in a N x 3 block (same as for atoms), used when [TrackedMode] = [Inline]. Rows must be ordered in the same way as in the [System%Atoms] block.
ModePath
Type: String Path to a .rkf file containing the modes which are to be tracked. Which modes will be refined is selected using the criteria from the [ModeSelect] block.)
ModeSelect
Type: Block Pick which modes to track from modes generated from Hessian or read from file.
FreqAndIRRange
Type: Float List cm-1 and km/mol True Specifies a combined frequency and IR intensity range within which all modes will be tracked. (First 2 numbers are the frequency range, last 2 numbers are the IR intensity range.)
FreqRange
Type: Float List cm-1 True Specifies a frequency range within which all modes will be tracked. (2 numbers: a upper and a lower bound.)
Full
Type: Bool False Track all modes.
HighFreq
Type: Integer Track the N modes with the highest frequencies.
HighIR
Type: Integer Track the N modes with the largest IR intensities.
IRRange
Type: Float List km/mol True Specifies an IR intensity range within which all modes will be tracked. (2 numbers: a upper and a lower bound.)
ImFreq
Type: Bool False Track all modes with imaginary frequencies.
LowFreq
Type: Integer Track the N modes with the lowest frequencies. (Includes imaginary modes which are recorded with negative frequencies.)
LowFreqNoIm
Type: Integer Track the N modes with the lowest non-negative frequencies. (Imaginary modes have negative frequencies and are thus omitted here.)
LowIR
Type: Integer Track the N modes with the smallest IR intensities.
ModeNumber
Type: Integer List Indices of the modes to track.
ScanModes
Type: Bool False If enabled an additional displacement will be performed along the new modes at the end of the calculation to obtain refined frequencies and IR intensities. Equivalent to running the output file of the mode tracking calculation through the AMS ModeScanning task.
ToleranceForBasis
Type: Float 0.0001 Convergence tolerance for the contribution of the newest basis vector to the tracked mode.
ToleranceForNorm
Type: Float 0.0005 Convergence tolerance for residual RMS value.
ToleranceForResidual
Type: Float 0.0005 Convergence tolerance for the maximum component of the residual vector.
TrackedMode
Type: Multiple Choice File [Inline, File, Hessian] Set how the initial guesses for the modes are supplied.
TrackingMethod
Type: Multiple Choice OverlapInitial [OverlapInitial, DifferenceInitial, FreqInitial, IRInitial, OverlapPrevious, DifferencePrevious, FreqPrevious, IRPrevious, HighestFreq, HighestIR, LowestFreq, LowestResidual] Set the tracking method that will be used.
UpdateMethod
Type: Multiple Choice JD [JD, D] Chooses the method for expanding the Krylov subspace: (D) Davidson or (JD) vdVorst-Sleijpen variant of Jacobi-Davidson.