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).

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. [11]

vdVorst-Sleijpen-Jacobi-Davidson:¶

This variant of the Jacobi-Davidson scheme from Sleijpen & vdVorst [11] 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. [11]

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.

Additional input parameters:¶

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
Description:	Input data for ModeTracking task.

Displacement

Type:	Float
Default value:	0.01
Description:	Step size (in Bohr) for finite difference calculation of frequencies and IR intensities during mode tracking iterations.

HessianGuess

Type:	Multiple Choice
Default value:	UFF
Options:	[Unit, File, UFF, Inline]
Description:	Sets how to obtain the guess for the Hessian used in the preconditioner.

HessianInline

Type:	Non-standard block
Description:	Initial guess for the (non-mass-weighted) Hessian in a 3N x 3N block, used when [HessianGuess] = [Inline].

HessianPath

Type:	String
Description:	Path to a .rkf file containing the initial guess for the Hessian, used when [HessianGuess] = [File].

MassWeighInlineMode

Type:	Bool
Default value:	True
Description:	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
Description:	Maximum number of allowed iterations.

ModeInline

Type:	Non-standard block
Recurring:	True
Description:	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
Description:	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
Description:	Pick which modes to track from modes generated from Hessian or read from file.

FreqAndIRRange

Type:	Float List
Unit:	cm-1 and km/mol
Recurring:	True
Description:	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
Unit:	cm-1
Recurring:	True
Description:	Specifies a frequency range within which all modes will be tracked. (2 numbers: a upper and a lower bound.)

Full

Type:	Bool
Default value:	False
Description:	Track all modes.

HighFreq

Type:	Integer
Description:	Track the N modes with the highest frequencies.

HighIR

Type:	Integer
Description:	Track the N modes with the largest IR intensities.

IRRange

Type:	Float List
Unit:	km/mol
Recurring:	True
Description:	Specifies an IR intensity range within which all modes will be tracked. (2 numbers: a upper and a lower bound.)

ImFreq

Type:	Bool
Default value:	False
Description:	Track all modes with imaginary frequencies.

LowFreq

Type:	Integer
Description:	Track the N modes with the lowest frequencies. (Includes imaginary modes which are recorded with negative frequencies.)

LowFreqNoIm

Type:	Integer
Description:	Track the N modes with the lowest non-negative frequencies. (Imaginary modes have negative frequencies and are thus omitted here.)

LowIR

Type:	Integer
Description:	Track the N modes with the smallest IR intensities.

ModeNumber

Type:	Integer List
Description:	Indices of the modes to track.

ScanModes

Type:	Bool
Default value:	False
Description:	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
Default value:	0.0001
Description:	Convergence tolerance for the contribution of the newest basis vector to the tracked mode.

ToleranceForNorm

Type:	Float
Default value:	0.0005
Description:	Convergence tolerance for residual RMS value.

ToleranceForResidual

Type:	Float
Default value:	0.0005
Description:	Convergence tolerance for the maximum component of the residual vector.

TrackedMode

Type:	Multiple Choice
Default value:	File
Options:	[Inline, File, Hessian]
Description:	Set how the initial guesses for the modes are supplied.

TrackingMethod

Type:	Multiple Choice
Default value:	OverlapInitial
Options:	[OverlapInitial, DifferenceInitial, FreqInitial, IRInitial, OverlapPrevious, DifferencePrevious, FreqPrevious, IRPrevious, HighestFreq, HighestIR, LowestFreq, LowestResidual]
Description:	Set the tracking method that will be used.

UpdateMethod

Type:	Multiple Choice
Default value:	JD
Options:	[JD, D]
Description:	Chooses the method for expanding the Krylov subspace: (D) Davidson or (JD) vdVorst-Sleijpen variant of Jacobi-Davidson.