Molecular Vibrations and IR/Raman Spectroscopy

A geometry optimization determines where the first derivative of the potential energy is small. A molecular vibration calculation examines the second derivative: it determines how the energy curves around that geometry, decomposes small nuclear motion into normal modes, and—when the calculator provides electronic response—assigns infrared intensities and Raman activities. These are different questions. Force convergence is necessary input to a vibration calculation; it is not a substitute for the curvature calculation.

The reader who needs the visible command path should use Calculate Molecular Vibrations. Script construction, validation, result fields, and executable sensitivity studies are kept in Molecular Vibrations with Tako Script. Those chapters apply the theory here without reintroducing it.

The local harmonic model

Let the Cartesian nuclear coordinates be collected in the $3N$-component vector $\mathbf{x}$, and let $\mathbf{x}_0$ denote the geometry about which the vibration is calculated. A Taylor expansion of the potential energy is

where $\mathbf{g}_0=\nabla E(\mathbf{x}_0)$ and

is the Cartesian Hessian. Indices $i,j$ label atoms and $\alpha,\beta$ label Cartesian directions. At an exact stationary point $\mathbf{g}_0=\mathbf{0}$, leaving a quadratic energy model. The harmonic approximation discards cubic and higher derivatives. It is consequently local: it describes sufficiently small oscillations in one basin, not bond dissociation, large-amplitude torsion, proton transfer, or barrier crossing.

Residual gradient matters even though the Hessian can be evaluated at any geometry. If $\mathbf{g}_0\ne\mathbf{0}$, the quadratic model is tilted. Low-frequency translations, rotations, and soft internal motions can mix, and a negative eigenvalue can describe continued relaxation rather than the curvature of a stationary species. A numerical frequency attached to a visibly nonstationary geometry is not physically meaningful.

Mass weighting and the normal-mode problem

Newton’s equations for a small displacement are

where $\mathbf{M}$ is diagonal with the atomic mass $m_i$ repeated for the three Cartesian components of atom $i$. Substituting a periodic displacement $\Delta\mathbf{x}=\mathbf{a}\exp(i\omega t)$ gives the generalized eigenproblem

Introducing mass-weighted coordinates $\mathbf{q}=\mathbf{M}^{1/2}\Delta\mathbf{x}$ converts this to an ordinary symmetric problem,

$\mathbf{F}$ is the mass-weighted Hessian, $\omega_k^2$ is its $k$th eigenvalue, and $\mathbf{e}_k$ is a mass-weighted eigenvector. With the convention $\mathbf{e}_k^{\mathrm T}\mathbf{e}_l=\delta_{kl}$, the corresponding Cartesian displacement is proportional to $\mathbf{M}^{-1/2}\mathbf{e}_k$. Light atoms therefore exhibit larger Cartesian amplitudes than heavy atoms for comparable mass-weighted components. The sign of an eigenvector is arbitrary: $\mathbf{e}_k$ and $-\mathbf{e}_k$ describe the same mode.

Frequency is reported as wavenumber,

usually in cm$^{-1}$. A positive eigenvalue gives a real $\omega_k$. A negative eigenvalue gives an imaginary angular frequency. Tako represents that case as a negative signed wavenumber in frequencies_cm_inv and also retains real/imaginary components in the structured result. The negative sign is a display convention for negative curvature, not a negative rate of oscillation.

Coordinate and unit conventions

The Cartesian Hessian depends on the coordinate unit used to take derivatives, while the physical frequencies do not after all conversions and masses are applied consistently. Tako accepts structure coordinates in ångström, converts the finite displacement to the engine’s internal Bohr convention, and returns frequencies in cm$^{-1}$. The raw eigenvectors are tied to the engine’s mass-weighting and normalization convention. Combining a Hessian, eigenvector, or force constant from another program with Tako arrays requires an explicit unit and normalization transformation; matching array dimensions is not sufficient.

Under a rigid rotation represented by an orthogonal Cartesian transformation $\mathbf{R}$, the unweighted Hessian transforms as $\mathbf{H}'=\mathbf{R}^{\mathrm T}\mathbf{H}\mathbf{R}$. Its individual Cartesian elements change, but the mass-weighted eigenvalues do not. This is why frequencies are rotationally invariant while printed Cartesian mode components depend on molecular orientation. A comparison based on one component such as “atom 3 moves in x” is fragile; compare the full displacement pattern in a common orientation or use internal-coordinate character.

For nondegenerate normalized modes, an overlap between calculations can be written

An overlap near one supports the same mass-weighted motion, independent of the arbitrary eigenvector sign. Before taking this dot product, atom order, isotope masses, coordinate frame, active-coordinate mask, and normalization must match. For a degenerate or nearly degenerate set, the individual eigenvectors can rotate freely within their subspace. Then compare the projector $\mathbf{P}=\sum_{k\in S}\mathbf{e}_k\mathbf{e}_k^{\mathrm T}$ or the singular values between the two mode subspaces rather than demanding a diagonal one-to-one overlap matrix.

External and internal modes

An isolated nonlinear molecule has three overall translations, three overall rotations, and $3N-6$ internal vibrations. A linear molecule has only two independent rotations and therefore $3N-5$ internal vibrations. In the ideal continuous problem, translation and rotation do not change the energy, so their Hessian eigenvalues are zero.

Tako diagonalizes the mass-weighted Cartesian Hessian without a separate translation/rotation projection step. The returned list therefore contains all $3N$ modes. For a well-converged nonlinear molecule, six should be near zero; for a linear molecule, five should be near zero. Finite force accuracy, finite displacement, incomplete geometry convergence, and floating-point algebra move these modes away from exactly zero and can give them either sign.

No universal wavenumber cutoff separates a harmless external-mode residue from a chemically meaningful imaginary mode. Classification uses several pieces of evidence together:

• whether the mode count and shape correspond to translation or rotation;
• whether the eigenvector is localized on an internal coordinate or reaction motion;
• whether tighter geometry convergence changes it;
• whether changing the finite-difference displacement changes it;
• whether the same mode persists with a method appropriate to the system.

For nearly linear molecules and floppy complexes, external and low-frequency internal motions can mix strongly. Treating every small negative number as a transition state or deleting every low mode as “rotation” are equally unreliable shortcuts.

Curvature and stationary-point classification

After external zero modes are identified, the sign pattern of the internal or vibrational subspace classifies local curvature:

A first-order saddle assignment also requires the imaginary eigenvector to follow the intended reaction coordinate. Counting one negative frequency without inspecting the displacement can misidentify a torsional instability, fragment rotation, or unrelated rearrangement as the desired transition state.

Conversely, a frequency calculation cannot establish that a local minimum is the global minimum, the dominant conformer, or thermodynamically populated. It classifies curvature only in the basin and electronic state supplied to the calculator.

Numerical Hessian from forces

Tako’s current GFN2 and MLIP paths do not expose analytic second derivatives to the vibration operation. The Hessian is constructed by displacing coordinates and differentiating forces. Because $F_{j\beta}=-\partial E/\partial x_{j\beta}$, a central two-point estimate is

where $\mathbf{u}_{i\alpha}$ is the unit displacement of atom $i$ along direction $\alpha$. For $N_a$ active atoms this stencil requires $6N_a$ displaced force evaluations, in addition to any reference evaluations used by response stages.

When the four-point option is selected, the derivative is

This uses $12N_a$ displaced force evaluations. Its truncation error falls more rapidly with $\delta$, but the additional points do not cure noisy forces, a nonstationary geometry, or an inadequate potential-energy model. The engine symmetrizes the assembled Hessian,

because exact mixed derivatives are symmetric while independently differenced matrix elements can disagree numerically.

Choosing the displacement

The displacement $\delta$ balances truncation and cancellation error. If it is too large, the force sample extends beyond the locally quadratic region and includes anharmonic curvature. If it is too small, the force difference in the numerator approaches calculator noise and floating-point resolution. Dividing that noisy difference by a smaller $\delta$ amplifies the error.

The current default is 0.01 Å. This is a starting scale, not proof of convergence. A useful numerical check repeats the calculation at a second physically reasonable displacement—often 0.005 or 0.02 Å—and compares the internal modes and their eigenvectors. Large systematic shifts indicate force noise, insufficient optimization, mode mixing, or anharmonicity. Close or nearly degenerate modes may rotate within their shared eigenspace even when the subspace itself is stable; compare the span or displacement patterns rather than requiring one-to-one eigenvector signs.

Higher-frequency localized stretches generally tolerate a smaller harmonic neighborhood than soft intermolecular modes. A single displacement may not be optimal for every coordinate, but Tako currently applies one Cartesian amplitude to the active set. Report it with the frequencies.

Constraints and the active Hessian

Tako’s current public vibration workflow selects active atoms, not arbitrary active Cartesian coordinates. An atom constrained in all three directions is excluded from the finite-displacement index set. Every remaining atom contributes all three displacement coordinates. Consequently, a fully fixed environment produces a partial Hessian for the active atoms, but a one- or two-axis Cartesian constraint does not remove only those axes from the Hessian. The atom remains in the displacement set and the constraint-adjusted forces affect the derivative matrix.

For a fully fixed environment, the omitted atoms cannot translate, rotate, relax, or couple dynamically to the active atoms. Frequencies can therefore shift relative to a full calculation, and the external-mode count no longer follows the simple $3N-6$ rule for the complete structure. Partial-axis constraints create a different constrained derivative problem and should not be described as a clean coordinate-deletion Hessian. Preserve the complete constraint mask and state which atoms, if any, were fixed in all directions. Freezing atoms merely to reduce cost is defensible only when the removed coupling is unimportant to the stated observable.

Zero-point energy

For a real harmonic normal mode, the quantum-mechanical energy levels are

The ground-state contribution summed over real modes is the zero-point energy,

Tako returns this quantity in Hartree and electronvolts. It is derived from the harmonic frequencies of the supplied model. It is not a complete thermal correction, entropy, free energy, or standard-state conversion. Imaginary modes do not define stable harmonic oscillators, and their absolute wavenumbers do not enter a minimum’s thermochemistry as real vibrations.

For a balanced reaction or barrier, a zero-point correction is a difference between consistently calculated states. For example,

The electronic method, state, geometry protocol, finite-difference settings, active coordinate policy, and treatment of low/imaginary modes must be compatible across the terms. The sign of the correction is not predetermined.

Isotope dependence

Within the Born–Oppenheimer harmonic model, isotope substitution changes $\mathbf{M}$ but not the electronic Hessian $\mathbf{H}$. Frequencies therefore change through mass weighting. For an ideal localized diatomic-like stretch with reduced mass $\mu$, $\omega\propto\mu^{-1/2}$. Real polyatomic modes are mixtures, so the shift depends on how strongly the substituted atom participates in the eigenvector.

Tako’s public structure format uses element identities and the engine’s associated masses; it does not currently expose a per-atom isotope-mass override in the vibration operation. Do not claim an isotope spectrum by relabeling peaks after the calculation. An isotope-specific prediction requires an implementation path that actually changes the mass matrix.

Infrared response

A harmonic frequency indicates an allowed nuclear motion, not whether the mode absorbs infrared radiation. In the electric-dipole approximation, a mode is IR active when the molecular dipole moment changes along its normal coordinate $Q_k$. The line strength is proportional to

Tako’s GFN2 vibration path performs a response calculation for requested infrared intensities and returns intensities_au, the magnitude of the static dipole, and the response calculation’s maximum-force diagnostic. The field name intensities_au is historical: the stored values follow the ASE (D/Å)² amu⁻¹ convention, while static_dipole_au is the scalar dipole norm in $e a_0$. The response is obtained from displaced electronic calculations and projected onto the normal modes. An MLIP supplies energy and forces but no electronic dipole response; when infrared output is requested with an MLIP, Tako logs that semi-empirical properties are required and returns the frequency result without an infrared object.

Symmetry provides selection rules. A centrosymmetric molecule has mutually exclusive ideal IR and Raman activity for gerade/ungerade modes, and a homonuclear diatomic has no dipole derivative. Numerical geometry asymmetry and finite differences can give a nominally forbidden mode a small residual intensity. Interpret weak lines against numerical convergence and symmetry rather than assigning physical meaning to every nonzero value.

The dipole derivative inherits two numerical layers: the displaced electronic response and the normal-mode transformation built from the force Hessian. An SCF state that changes discontinuously across $+\delta$ and $-\delta$, an incorrect charge/spin assignment, or a mode that rotates within a near-degenerate pair can perturb the intensity even when the frequency appears stable. Intensity convergence should therefore test the electronic settings and displacement, not only the final broadened plot.

Absolute experimental absorbance also depends on concentration, path length, population, line shape, and instrument. Tako’s stored IR intensities are calculated response strengths in the ASE convention; they are not a Beer–Lambert absorbance trace. Normalizing a plot to its strongest peak discards absolute scale and must be stated when comparisons are made.

Raman response

Raman scattering depends on the derivative of the electronic polarizability tensor $\boldsymbol\alpha$ along a normal coordinate. In the nonresonant Placzek approximation, define the isotropic derivative

and anisotropic combinations collected in $\gamma_k'^2$. A commonly used orientation-averaged activity has the form

Tako’s GFN2 Raman path evaluates the electronic response and returns activities, absolute-intensity data, and invariants containing alpha2, gamma2, and delta2. These are method- and convention-dependent response quantities. Raman activity is not identical to a measured peak height: experimental intensity also depends on excitation frequency, temperature, polarization geometry, instrumental response, and line broadening.

An MLIP does not provide a polarizability derivative. Requested Raman output is therefore skipped with an explanatory log message while frequencies remain available.

The nonresonant Placzek model assumes that the laser is sufficiently far from an electronic transition for a static or slowly varying polarizability derivative to describe scattering. Resonance Raman enhancement, excited-state vibronic coupling, and laser-frequency dependence are outside this approximation. In the current runtime, both activities and absolute_intensities are the same $45\alpha^2+7\gamma^2$ activity array in Å$^4$/amu because the antisymmetric-anisotropy weight is zero. The latter name does not mean that laser frequency, temperature, population, or experimental collection geometry has been applied. For comparison with experiment, retain the invariants and state any later intensity model explicitly.

Polarization selection rules carry additional structural information. Isotropic and anisotropic invariants distinguish totally symmetric character from depolarized scattering contributions. A powder or solution average is not the same observable as a fixed single-molecule orientation. Tako returns invariants but does not infer the experimental polarization geometry from a structure file.

Stick data and plotted spectra

The fundamental calculation produces discrete mode frequencies and response strengths. A continuous-looking spectrum is a visualization obtained by broadening those sticks with a chosen line shape and width. Broadening does not improve the electronic or vibrational calculation and does not create resolving power absent from the underlying modes.

When comparing with experiment, report whether the comparison uses raw harmonic sticks or a broadened plot, the broadening kernel and width, any frequency scaling, temperature, and normalization. Tako’s returned vibration/IR/Raman objects contain the calculated mode data; a viewer presentation is not additional computed physics.

Harmonic frequencies are commonly shifted from experimental fundamentals by electronic-structure error and anharmonicity. An empirical scale factor is method-, basis-, and observable-specific. Tako does not apply a frequency-scaling keyword or a universal correction in the vibration result. If a scale factor is used in analysis, retain both raw and transformed frequencies and cite the factor’s validation source.

Tako execution sequence

The vibration operation first constructs frequencies and modes. It emits a cumulative partial result when those modes are ready. If infrared is requested and supported, a separate GFN2 response pass follows and another partial result is emitted. Raman is processed afterward in its own response pass. The final promise resolves only after all requested supported stages finish.

The sequence matters for failure interpretation. A late infrared or Raman error can occur after a valid frequency partial result exists, but a partial workspace artifact is not the same as a successfully resolved all-requested-properties result. Cancellation exposes no restartable finite-difference checkpoint. A later call begins a new operation and repeats its required displaced evaluations.

GFN2 is the current route to frequencies plus IR/Raman electronic response. MLIPs can produce the force-derived Hessian and frequencies but skip electronic response objects. Calculator capability, numerical completion, stationary-point quality, and spectroscopic accuracy are separate conditions.

Computational cost

For $N_a$ active atoms, the two-point Hessian requires approximately $6N_a$ force evaluations and the four-point Hessian approximately $12N_a$. Each requested electronic response channel can add its own displaced evaluations, so frequencies plus IR plus Raman is substantially more expensive than frequencies alone. The actual cost is

For browser-scale systems the calculator evaluations normally dominate until the dense $3N_a\times3N_a$ Hessian and eigensolver become large. Reducing the active set with constraints also changes the physical Hessian, so it is not a free performance optimization. Request only the response channels required by the scientific question.

Numerical validation

A defensible vibration result is supported by a convergence matrix rather than one run:

Agreement in frequency alone can hide a swapped mode assignment. Compare eigenvector patterns, localization, symmetry character, and ordering. Near-degenerate modes can mix under tiny perturbations; compare their combined subspace and total response strength rather than insisting on a fixed individual labeling.

Separating error sources

Four error classes can produce a shifted or imaginary mode, and they require different remedies:

1. Geometry error. The reference structure is not stationary under the same calculator. Tighten or repair the optimization and preserve its state settings.
2. Differentiation error. The displacement or stencil does not resolve the force derivative. Test $\delta$ and the two-/four-point choice.
3. Calculator error or noise. The potential gives inaccurate curvature or forces fluctuate at the required scale. Validate against a matched higher-level or experimental reference; more stencil points cannot repair the model.
4. Harmonic-model error. The local quadratic approximation is physically inadequate. Large-amplitude or anharmonic treatment is required rather than tighter numerical differentiation.

Changing several controls simultaneously can make a frequency look better without identifying which error was reduced. A convergence study changes one numerical element at a time, keeps the geometry/electronic state traceable, and compares both frequencies and mode character.

Soft modes and conformational coordinates

Low-frequency internal modes contribute disproportionately to entropy in a rigid-rotor/harmonic-oscillator treatment, yet they are also the modes least well described as quadratic oscillators. Torsions, hindered rotations, fragment librations, and weak intermolecular translations can span multiple minima. Tako reports their harmonic curvature but does not replace them with hindered-rotor partition functions or quasi-harmonic corrections.

This boundary matters even when the immediate purpose is minimum classification. A small positive soft mode supports positive local curvature, but its numerical wavenumber may be very sensitive to conformation, solvent environment, constraints, and displacement. Stable stiff stretching modes alone do not make a harmonic ZPE or later thermochemical correction high-precision.

Relation to thermochemistry

Frequencies are inputs to common rigid-rotor/harmonic-oscillator thermal corrections, but Tako’s vibration result currently returns the harmonic ZPE rather than a complete temperature- and pressure-dependent thermochemistry table. Translational, rotational, vibrational, and electronic partition functions require molecular mass, moments of inertia, symmetry number, electronic degeneracy, temperature, standard state, and a policy for low or imaginary modes.

Standard thermochemistry treatments develop these terms and their program-specific controls. Tako’s vibration result does not evaluate them. If thermal enthalpy, entropy, or free energy is calculated downstream, document the equations and conventions separately and preserve the raw frequencies. Standard-state conversions can be comparable to the chemical energy difference of interest, especially for association reactions.

For transition-state theory, one imaginary reaction-coordinate mode is excluded from the stable harmonic partition function and treated through the rate model. Simply taking its absolute frequency and adding it to ZPE or entropy produces a different, generally invalid, calculation.

Reading the result

The base vibrations object contains signed frequencies, real/imaginary frequency components, harmonic energies, zero-point energy, and Cartesian mode arrays obtained as $\mathbf{M}^{-1/2}\mathbf e_k$ from normalized mass-weighted eigenvectors. The optional infrared and raman objects each retain their own vibration payload alongside response quantities. This duplication records the modes used in each response pass and should not be flattened without checking consistency.

Mode arrays are numerical displacement patterns, not trajectories at a physical temperature or time step. Their arbitrary sign and normalization mean that animation amplitude is a viewer choice. Use the relative atomic directions to assign the motion.

Acceptance of a molecular minimum normally requires:

• a converged input geometry under a recorded force convention;
• the expected near-zero external modes;
• no persistent meaningful imaginary internal mode;
• stable important frequencies under displacement/stencil checks;
• a calculator and electronic state appropriate to the molecule;
• explicit separation of raw sticks, response strengths, and any plotted transformation.

Acceptance of a transition-state candidate instead requires exactly one meaningful imaginary internal mode aligned with the intended reaction coordinate, with other internal modes real within the numerical protocol.

Limitations of the harmonic model

The harmonic approximation treats each normal coordinate as an independent quadratic oscillator. Real molecular vibrations exhibit anharmonic level spacing, mode coupling, resonances, tunneling, large-amplitude internal rotation, and temperature-dependent conformational sampling. These effects are most important for X–H stretches, floppy intermolecular modes, proton transfer, and near-degenerate vibrations.

Tako’s current vibration operation does not compute anharmonic force fields, VPT2 corrections, rovibrational structure, resonance treatment, temperature-dependent populations, or explicit-solvent vibrational ensembles. Implicit solvation changes the electronic model but does not replace sampling of solvent configurations. A harmonic spectrum can support assignments and comparative trends within its domain; it should not be presented as a complete reproduction of an experimental spectrum without stating these omissions.

Reporting

Report the molecular structure and electronic state, calculator/model and asset provenance, dispersion and solvation, geometry-optimization threshold and constraint mask, finite-difference displacement and stencil, active atom count, signed frequencies, treatment of external and imaginary modes, ZPE convention, requested response channels, and whether any scaling or broadening was applied after calculation. For a key band, include a mode description or image rather than only its index, because indices can change when close modes reorder.

References

1. E. B. Wilson, J. C. Decius, P. C. Cross, Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra, McGraw-Hill, New York (1955).
2. D. A. Long, The Raman Effect: A Unified Treatment of the Theory of Raman Scattering by Molecules, Wiley, Chichester (2002).
3. P. Atkins, R. Friedman, Molecular Quantum Mechanics, Oxford University Press, Oxford.
4. C. Bannwarth, S. Ehlert, S. Grimme, “GFN2-xTB—An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions,” J. Chem. Theory Comput. 15, 1652–1671 (2019).
5. A. H. Larsen et al., “The Atomic Simulation Environment—A Python Library for Working with Atoms,” J. Phys.: Condens. Matter 29, 273002 (2017).