nitrogen.vpt

Vibrational perturbation theory and harmonic oscillator methods

nitrogen.vpt.a_ap_matrix(n)

Calculate the matrix representations of the annihilation and creation operators, \(a\) and \(a^\dagger\).

Parameters

n (int) – The number of states. The last state is \(\vert n-1 \rangle\).

Returns

• a (ndarray) – The (n, n) matrix representation of \(a\).

• ap (ndarray) – The (n, n) matrix representation of \(a^\dagger\).

nitrogen.vpt.analyze_vertff(w, f)

Analyze a vertical harmonic force field

Parameters

• w ((n,) array_like) – The initial state harmonic frequencies

• f ((...,) array_like) – The vertical derivative array up to at least quadratic

Returns

• omega (ndarray) – The absolute magnitude of the upper state frequencies

• sigma (ndarray) – The imaginary coefficient of the upper state frequencies. (1 or minus i)

• d (ndarray) – The location of the upper state stationary point with respect to the lower state coordinates

• v0 (float) – The stationary point energy.

nitrogen.vpt.autocorr_cubic_initial_firstorder(V, f, n, t)

Calculate the first-order corrections to the quadratic auto-correlation function due to cubic force constants in the initial state.

Parameters

• V (ndarray) – The derivative array of the initial state, including up to cubic derivatives.

• f (ndarray) – The derivative array of the final state, including up to quadratic derivatives.

• n (int) – The number of coordinates

• t (ndarray) – The time array.

Returns

C1 – The first-order correction relative to C0.

Return type

ndarray

Notes

The returned value C1 is the ratio of the first-order correction to the zeroth-order quadratic correlator. The total result is \(\left[ 1+C_1(t)\right] \times C_0(t)\).

nitrogen.vpt.autocorr_cubic_initial_firstorder_linearHT(V, f, mu, n, t)

Calculate the first-order corrections to the quadratic auto-correlation function due to cubic force constants in the initial state and include linear Herzberg-Teller terms.

Parameters

• V (ndarray) – The derivative array of the initial state, including up to cubic derivatives.

• f (ndarray) – The derivative array of the final state, including up to quadratic derivatives.

• mu (array_like) – The dipole surface, including up to gradients.

• n (int) – The number of coordinates

• t (ndarray) – The time array.

Returns

C1LHT – The dipole autocorrelation function relative to C0.

Return type

ndarray

Notes

nitrogen.vpt.autocorr_linear(w, f, t)

Calculate the vacuum state autocorrelation function for propagation on a linear potential energy surface.

Parameters

• w (array_like) – The harmonic frequency (in energy units) of each mode.

• f (array_like) – The derivative array, including at least first derivatives.

• t (array_like) – The time array, in units where \(\hbar = 1\). (Alternatively, the t array can be identified with \(t/\hbar\).)

Returns

C – The autocorrelation function, \(C(t)\).

Return type

ndarray

See also

spech_fft

Calculate the spectrum of an autocorrelation function

Notes

A linear potential is separable into 1-D components, so the total autocorrelation function is a product with a factor for each mode equal to

\[C_i(t) = \frac{1}{(1 + i \omega t )^{1/2}} \exp[-(f^2 t^2/24)(6 + i \omega t)]\]

times a final phase factor of \(\exp[-i f_0 t]\), where \(f_0\) is the energy offset equal to f[0].

nitrogen.vpt.autocorr_linearHT(w, f, mu, t)

Calculate the harmonic vacuum autocorrelation function including linear Herzberg-Teller dipole terms.

Parameters

• w (array_like) – The harmonic frequency (in energy units) of each mode.

• f (array_like) – The derivative array of the final state, including at least second derivatives.

• mu (array_like) – The derivative array of the dipole function, including at least gradients.

• t (array_like) – The time array.

Returns

CHT – The dipole autocorrelation function relative to C0.

Return type

ndarray

Notes

The returned value CHT is the ratio of the dipole correlation function relative to the quadratic correlator, C0. The total result is \(C_{HT}(t)\times C_0(t)\).

See also

autocorr_quad

Quadratic autocorrelation function.

nitrogen.vpt.autocorr_quad(W, f, t, method='direct')

Calculate the vacuum state autocorrelation function for propagation on a quadratic potential energy surface.

Parameters

• W ((n,) or (n,n) array_like) – The harmonic frequency (in energy units) of each mode, or the inverse mass tensor.

• f (array_like) – The derivative array, including up to at least second derivatives.

• t (array_like) – The time array, in units where \(\hbar = 1\). (Alternatively, the t array can be identified with \(t/\hbar\).)

• method ({'direct','integral','integral_log'}) – The calculation method. See Notes

Returns

C – The autocorrelation function, \(C(t)\).

Return type

ndarray

See also

corr_quad_recursion_elements

Calculate quadratic correlator recursion coefficients

spech_fft

Calculate the spectrum of an autocorrelation function

Notes

For method = ‘direct’, a direct expression based on a discontinuity-free BCH disentangling formula is used.

For method = ‘integral’, an alternative method is used to first calculate the logarithmic derivative of \(C(t)\). This is numerically integrated by a cumulative version of Simpson’s rule and then exponentiated.

For method = ‘integral_log’, the integrated logarithm is returned directly, without exponentiation. That is, the branch-cut discontinuity-free logarithm of \(C(t)\) is returned.

For the integral methods, a sufficiently small time-step in the t array is required for accurate results. The direct method does not rely in numerical integration.

nitrogen.vpt.autocorr_quad_finiteT(w, f, t, beta, method='direct')

Calculate the thermal correlation function for propagation on a quadratic potential energy surface.

Parameters

• w (array_like) – The harmonic frequency (in energy units) of each mode.

• f (array_like) – The derivative array, including up to at least second derivatives.

• t (array_like) – The time array, in units where \(\hbar = 1\). (Alternatively, the t array can be identified with \(t/\hbar\).)

• beta (float) – The value of \(\beta = 1/kT\) in inverse energy units.

• method ({'direct'}) – The calculation method. See Notes.

Returns

C – The thermal autocorrelation function, \(C(t)\).

Return type

ndarray

See also

autocorr_quad

Calculate the zero-temperature correlation function

spech_fft

Calculate the spectrum of an autocorrelation function

Notes

For method = ‘direct’ (currently the only method), a direct evaluation of the thermal trace

\[C(t, \beta) = \frac{1}{Z_0(\beta)} \text{Tr}\left[ e^{(+it-\beta)H_0} e^{-it H_1} \right]\]

is performed using a stable, discontinuity-free method.

nitrogen.vpt.calc_rectilinear_modes(hes, mass, hbar=None, norm='dimensionless')

Calculate the rectilinear normal modes and energies

Parameters

• hes (array_like) – The (3*N,3*N) Cartesian Hessian matrix.

• mass (array_like) – The N masses

• hbar (float, optional) – The value of \(\hbar\). If None (default), NITROGEN units will be assumed.

• norm ({'dimensionless', 'mass-weighted','both'}) – The normalization convention of the displacement vectors.

Returns

• w ((N,) ndarray) – The harmonic frequencies, in energy units.

• R ((3*N,3*N) ndarray) – Each column of R is the displacement vector for the corresponding normal mode. If norm is ‘’both’’, then T and L are returned in that order.

Notes

For norm = ‘dimensionless’, the displacement vectors are those of the non-mass-weighted Cartesian coordinates with respect to dimensionless, normalized coordinates. In these coordinates, the potential energy surface is

\[V = \sum_i \frac{1}{2} \omega_i q_i^2\]

For norm = ‘mass-weighted’, the displacement vectors are those for the mass-weighted Cartesian coordinates and equal the eigenvectors of the mass-weighted Hessian, i.e. the traditional \(\mathbf{L}\) array.

nitrogen.vpt.corr_quad_ratio_table(r, S, T, nmax)

Calculate correlation function ratios using the quadratic recursion elements.

Parameters

• r (ndarray) – Quadratic recursion coefficients.

• S (ndarray) – Quadratic recursion coefficients.

• T (ndarray) – Quadratic recursion coefficients.

• nmax (int) – The maximum total quantum number.

Returns

• Imn ((…,ns,ns) ndarray) – The correlator ratios. ns is the number of states.

• qns ((ns, n) ndarray) – The quantum number table for n modes.

See also

corr_quad_recursion_elements

Calculate the recursion coefficients.

corr_quad_ratio_table_edge

Calculate the edge of the ratio table.

nitrogen.vpt.corr_quad_ratio_table_edge(r, S, T, nmax)

Calculate correlation function ratios of type <m|…|0> using the quadratic recursion elements.

Parameters

• r (ndarray) – Quadratic recursion coefficients.

• S (ndarray) – Quadratic recursion coefficients.

• T (ndarray) – Quadratic recursion coefficients.

• nmax (int) – The maximum total quantum number.

Returns

• Im0 ((…,ns) ndarray) – The edge correlator ratios. ns is the number of states.

• qns ((ns, n) ndarray) – The quantum number table for n modes.

See also

corr_quad_recursion_elements

Calculate the recursion coefficients.

corr_quad_ratio_table

Calculate all correlator ratios.

nitrogen.vpt.corr_quad_ratio_table_rectangular(r, S, T, mmax, nmax)

Calculate correlation function ratios using the quadratic recursion elements for a rectangularly shaped table.

Parameters

• r (ndarray) – Quadratic recursion coefficients.

• S (ndarray) – Quadratic recursion coefficients.

• T (ndarray) – Quadratic recursion coefficients.

• mmax (int) – The maximum total quantum number (left index).

• nmax (int) – The maximum total quantum number (right index).

Returns

• Imn ((…,nsm,nsn) ndarray) – The correlator ratios.

• qns ((ns, n) ndarray) – The quantum number table for n modes. ns is the greater of nsm and nsn.

See also

corr_quad_recursion_elements

Calculate the recursion coefficients.

corr_quad_ratio_table_edge

Calculate the edge of the ratio table.

corr_quad_ratio_table

Calculate a square ratio table.

nitrogen.vpt.corr_quad_recursion_elements(W, f, t)

Calculate the correlation function recursion coefficients for a quadratic Hamiltonian.

Parameters

• W ((n,) or (n,n) array_like) – If a 1-d vector, the harmonic frequencies (in energy units) defining the dimensionless normal coordinates. If a 2-d array, then the effective inverse mass tensor.

• f (array_like) – The derivative array of the propagating surface. This must have at least second derivatives. (Note the factorial scaling rules of derivative arrays.)

• t (array_like) – The scaled time array, \(t/\hbar\).

Returns

r, S, T – The coefficients arrays

Return type

ndarray

Notes

Following Ref. [FT1989] , the correlation functions obey the recursion relations given by

\[ \begin{align}\begin{aligned}\langle m_i + 1 \vert \vert \cdots \rangle &= \frac{1}{\sqrt{m_i+1}} \left[ -r_i \langle \cdots \vert \vert \cdots \rangle + \sum_j S_{ij} \sqrt{n_j} \langle \cdots \vert \vert n_j - 1 \rangle - T_{ij} \sqrt{m_j} \langle m_j - 1 \vert \vert \cdots \rangle \right]\\\langle \cdots \vert \vert n_i + 1 \rangle &= \frac{1}{\sqrt{n_i+1}} \left[ -r_i \langle \cdots \vert\vert\cdots\rangle + \sum_j -T_{ij} \sqrt{n_j} \langle \cdots \vert \vert n_j -1 \rangle + S_{ij} \sqrt{m_j} \langle m_j - 1 \vert \vert \cdots \rangle \right]\end{aligned}\end{align} \]

where some simplications from the most general expressions have been applied given that the propagation is unitary. The matrices \(S\) and \(T\) are symmetric.

References

FT1989

F. M. Fernandez and R. H. Tipping, “Multidimensional harmonic oscillator matrix elements”. J. Chem. Phys., 91, 5505 (1989). https://doi.org/10.1063/1.457553

nitrogen.vpt.cubic_firstorder_vacuum(V, n)

Calculate the first-order coefficients of the vacuum ground state from cubic anharmonic force constants.

Parameters

• V (ndarray) – The scaled derivative array.

• n (int) – The number of coordinates.

Returns

c – The first order coefficients.

Return type

ndarray

Notes

The coefficient array c has the same format and sorting as derivative arrays, where the quantum numbers replace the deriative multi-index.

nitrogen.vpt.cubic_gradient_estimate(Vi, Vf, n)

Calculate the estimated spectral shift due to first-order cubic corrections of the lower state vacuum with the upper state gradient.

Parameters

• Vi (ndarray) – The derivative array of the lower state containing at least cubic derivatives.

• Vf (ndarray) – The deriative array of the upper state containing at least gradients.

• n (int) – The number of coordinates.

Returns

deltaE – The approximate spectral shift.

Return type

float

Notes

The shift is estimated as the product of the first-order displacement of the initial wavefunction with the final state gradient,

\[\Delta E \approx \sum_i \langle q_i \rangle f_i,\]

where \(f_i\) is the final state gradient. The coordinate displacements are evaluated to first-order as

\[\langle q_i \rangle \approx -\frac{1}{4 \omega_i} \left[\phi_{iii} + \sum_{j \neq i} \phi_{ijj} \right],\]

where \(\phi_{ijk}\) are the unscaled cubic derivatives.

nitrogen.vpt.linear_on_general(L, C, n)

Calculate the action of a general linear coordinate polynomial on a general HO expansion.

Parameters

• L (ndarray) – The derivative array, including at least gradients

• C (ndarray) – The coefficient array

• n (int) – The number of coordinates

Returns

Z – The coefficient array of L|C>

Return type

ndarray

nitrogen.vpt.p_matrix(n, m=1)

Calculate the matrix representation of \(p^m\), where \(p = -i \sqrt{\frac{1}{2}}(a - a^\dagger)\).

Parameters

• n (int) – The number of states. The last state is \(\vert n-1 \rangle\).

• power (int, optional) – The power of \(p\). The default is 1. Must be a non-negative integer.

Returns

p – The (n,`n`) matrix representation of \(p^m\).

Return type

ndarray

nitrogen.vpt.q_matrix(n, m=1)

Calculate the matrix representation of \(q^m\), where \(q = \sqrt{\frac{1}{2}}(a + a^\dagger)\).

Parameters

• n (int) – The number of states. The last state is \(\vert n-1 \rangle\).

• power (int, optional) – The power of \(q\). The default is 1. Must be a non-negative integer.

Returns

q – The (n,`n`) matrix representation of \(q^m\).

Return type

ndarray

nitrogen.vpt.cfourvib

CFOUR vibrational file processing and interface routines.

nitrogen.vpt.cfourvib.QUADRATURE2xyz(filename, elements, out='out.xyz', comment='')

Generate a .xyz file from a CFOUR QUADRATURE file.

Parameters

• filename (str) – The QUADRATURE file path.

• elements (array_like) – The element labels, e.g. [‘H’,’O’,’H’]

• out (str, optional) – The output file. The default is “out.xyz”.

• comment (str, optional) – The .xyz file comment string. The default is None.

Returns

Return type

None.

Notes

Note that the standard CFOUR QUADRATURE file uses bohr units, and the default .xyz file uses Angstroms units. This unit conversion is performed in this function.

nitrogen.vpt.cfourvib.read_QUADRATURE(filename, use_bohr=False, is_linear=False)

Parse a CFOUR QUADRATURE file.

Parameters

• filename (str) – The QUADRATURE file path.

• use_bohr (bool, optional) – If True, return displacements and geometry in bohrs. The default is False.

• is_linear (bool, optional) – If True, 3*N-5 normal modes will be parsed instead of 3*N-6. The default is False.

Returns

• freq ((nvib,) ndarray) – The harmonic frequencies

• T ((3*natom,nvib) ndarray) – The Cartesian displacement vectors, in Angstroms, of each dimensionless normal mode.

• ref_geo ((3*natom,) ndarray) – The reference Cartesian geometry in Angstroms.

nitrogen.vpt.cfourvib.read_cubic(filename, nvib=None, offset=7)

Read a CFOUR cubic force constants text file

Parameters

• filename (str) – The file path

• nvib (int, optional) – The number of vibrational modes. If None, the total will be inferred from the input file.

• offset (int, optional) – The mode numbering offset from zero. The default is 7, which is usually what CFOUR files require (6 for rot-trans modes and 1 for zero-indexing).

Returns

F – The scaled derivative including up to cubic derivatives. (The zeroth, first, and second derivatives are zeroed.)

Return type

ndarray

Notes

The deriative array is returned in standard scaled format, i.e. the derivatives are divided by the factorial of their respective multi-index.

nitrogen.vpt.cfourvib.sample_QUADRATURE(filename, num_sample, use_bohr=False, rng_seed=None, stddev=1.0, is_linear=False)

Generate sample points using a CFOUR QUADRATURE file.

Parameters

• filename (str) – The QUADRATURE file path.

• num_sample (integer) – The number of sample points.

• use_bohr (bool, optional) – If True, return Cartesian coordintes in bohrs. The default is False.

• rng_seed (integer, optional) – The RNG seed.

• stddev (float, optional) – The standard deviation of the sampled dimensionless normal coordinates, \(q\).

• is_linear (bool, optional) – Set to True if linear (i.e. 3*N-5 normal modes in the QUADRATURE file). Default is False.

Returns

X – The sampled geometries

Return type

(3*N, num_sample) ndarray

Notes

The harmonic oscillator ground-state wavefunction in terms of \(q\) is \(\psi(q) \propto \exp[-q^2/2]\), and the ground-state density is \(\vert \psi \vert ^2 \propto \exp[-q^2]\). These distributions have standard deviations of \(1\) and \(\sqrt{\frac{1}{2}} \approx 0.707\), respectively. Therefore, setting the option stddev = 0.707 will sample the harmonic ground-state density exactly. The default option (stddev = 1.0) samples the ground-state amplitude, which is \(\sqrt{2}\) wider.