Derivatives w.r.t. geometrical stuffs (qcip_tools.derivatives_g)

Tools to help the manipulation of derivatives of the energy with respect to cartesian (G) or normal mode (N) coordinates.

Theory

The displacement around the equilibrium position, \(\mathbf{x}_i=\mathbf{x}'_i-\mathbf{x}_{i,0}\), may be expressed as a Taylor series to describe the potential energy surface:

\[V(\mathbf{x}) = V_0 + \sum_i F_i\mathbf{x}_i + \frac{1}{2!}\,\sum_{ij}H_{ij}\mathbf{x}_i\mathbf{x}_j + \frac{1}{3!}\,\sum_{ij}F_{ijk}\mathbf{x}_i\mathbf{x}_j\mathbf{x}_k+ \ldots\]

where \(\mathbf{x}\in\mathbf{R}^{3N}\), \(F_i=\left.\frac{\partial V}{\partial\mathbf{x}_i}\right|_{\mathbf{x}=0}\) is the gradient (so the forces, wich are close to 0 if the system is in equilibrium), \(H_{ij}=\left.\frac{\partial^2 V}{\partial\mathbf{x}_i\mathbf{x}_j}\right|_{\mathbf{x}=0}\) is the hessian (force constant matrix) and \(F_{ijk}\) is the cubic force constants matrix (and so on). The sums runs over the \(3N\) degrees of freedom of \(N\) atoms.

Normal modes

By settings the Lagrange equations of the harmonic oscillator,

\[\begin{split}\begin{align} &T(\mathbf{x}) = \frac{1}{2} \sum_{ij} m_{ij}\,\dot{x}_i\dot{x}_j, \\ &V(\mathbf{x}) = \frac{1}{2}\sum_{ij} H_{ij}\,x_i\,x_j, \end{align}\end{split}\]

and \(L = T-V\). Taking advantage of the fact that the Hessian and the mass matrices are diagonals, \(H_{ij}=H_{ii}\delta_{ij}\) and \(m_{ij}=m_i\delta_{ij}\), the \(3N\) Lagrange equation are written

\[\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{x}_i}\right) - \left(\frac{\partial L}{\partial x_i}\right) = \sum_j m_{ij}\,\overset{\cdot\cdot}{x}_j + H_{ij}\,x_{j} = 0,\]

or, in a matrix form,

\[\mathbf{m}\,\overset{\cdot\cdot}{\mathbf{X}} + \mathbf{H}\,\mathbf{X} = 0.\]

One can assume, in the case of small displacements, a linear harmonic solution, \(\mathbf{X}=\mathbf{X}_0\,\cos{(\mathbf{\omega}\,t)}\), so that the previous equation reduce to

\[\mathbf{H}\,\mathbf{X}_0 = \mathbf{m}\,\mathbf{\omega}^2\,\mathbf{X}_0.\]

This is reduced to an eigenvalue problem by setting mass-weighted coordinates, \(\mathbf{Q}=\mathbf{m}^{\frac{1}{2}}\,\mathbf{X}_0\), which gives

\[(\mathbf{m}^{-\frac{1}{2}}\mathbf{H}\,\mathbf{m}^{-\frac{1}{2}})\,\mathbf{Q} = \mathbf{\omega}^2\,\mathbf{Q},\]

or,

(1)\[\mathbf{Q}^\dagger\,\mathbf{H}^m\,\mathbf{Q} = \omega^2\]

where \(\mathbf{H}^m=\mathbf{m}^{-\frac{1}{2}}\mathbf{H}\,\mathbf{m}^{-\frac{1}{2}}\) is the mass-weighted Hessian. The \(\mathbf{Q}\)’s are therefore eigenvectors of the mass-weigted hessian, also called normal modes (which form a complete and orthogonal basis), while the \(\mathbf{\omega}\)’s are the eigenvalues (vibrational frequencies). There is 6 (or 5, if the system is linear) zeros eigenvalues for the translation and rotation modes, while the other normal modes describe the vibrations.

Computation of normal modes

The computation takes place in the MassWeightedHessian class.

Warning

Modes are not decontamined from translations and rotations.

  1. Compute the mass weighted Hessian, \(\mathbf{H}^m\), as

    \[H^m_{i\alpha,j\beta} = \frac{1}{\sqrt{m_i\,m_j}}\,H_{i\alpha,j\beta},\]

where \(i\) and \(j\) are the number of the atom, while \(\alpha\) and \(\beta\) ar the cartesian coordinates of the atoms. Masses, \(m_i\) and \(m_j\), are in AMU.

  1. Compute the \(3N\) eigenvalues (\(\omega^2\)) and egeinvectors (\(\mathbf{Q}\)) of \(\mathbf{H}^m\). The explicit set of equations to solve, obtained from (1), is written

    (2)\[\omega_p^2\,\delta_{pq} = \sum_{i\alpha,j\beta} Q_{i\alpha,p}\,\frac{H_{i\alpha,j\beta}}{\sqrt{m_i\,m_j}}\,Q_{j\beta,q}.\label{eq:v1}\]

Order them from the lowest to the largest eigenvalue. Eventually discard the 6 (or 5) first lowest ones.

  1. Compute frequencies:

    \[\nu_p = \sqrt{\frac{\omega^2_p}{\mu}},\]

where \(\mu\) is the conversion factor from AMU to electron mass (about 1822, see here) and \(p\) is the normal mode. With the conversion, \(\nu_p\) is therefore in atomic units.

  1. Compute the carthesian displacements. From (2), they are defined by

    \[D_{p,i\alpha} = \frac{1}{\sqrt{m_i\,\mu}} Q_{i\alpha,p}.\]

With the conversion, displacement are therefore in \(\text{Å}\,m_e^{-1/2}\), where \(m_e\) is the electron mass (so atomic unit of mass).

Note

The normal modes (normal_modes) and cartesian displacements (displacements) are stored in the form \(p,i\alpha\).

Computation of thermochemistry data

See this excellent white paper from Gaussian on thermochemistry, which explains step by step how the thermochemistry quantities are computed.

Since the vibrational frequencies are needed, the MassWeightedHessian object also provide the methods for thermochemistry. Note that currently, the heat capacity is not implemented.

Note

The ZPVA contribution can be computed separately with compute_zpva(). Therefore it is not included in the computation of internal energy, enthalpy of free Gibbs energy.

API documentation

class qcip_tools.derivatives_g.BaseGeometricalDerivativeTensor(spacial_dof, trans_plus_rot=0, representation='', components=None)

Base for all derivatives of the energy with respect to geometrical parameters

to_name()

Get the name of the tensor from its representation, by using REPRESENTATION.

Return type:

str

class qcip_tools.derivatives_g.MassWeightedHessian(molecule, cartesian_hessian=None, scale=1.0)

To be exact, this is not a derivatives of the energy. So it is not a child of derivatives.Tensor !

Parameters:
ENERGY_IN_AU_CONVERSION = 3.808798847133371e-07

Convert energy from J/mol to Hartree (per particle)

HARTREE_TO_WAVENUMBER_CONVERSION = 219474.6313632

Convert Hartree to wavenumber (cm⁻¹)

INERTIA_CONVERSION = 1.6605390666000002e-47

Convert inertia from UMA*A² to kg*m²

MASS_CONVERSION = 1.6605390666e-27

Convert mass from AMU to kg

VIB_ENERGY_CONVERSION = 6579683920501762.0

Convert hartree to hertz

compute_enthalpy(temperature=298.15)

Compute the value of the different contribution (translation, rotation, vibration) to the enthalpy, IN HARTREE !

Note

Since \(H=U+RT\), the code just adds \(RT\) to the translational contribution to internal energy, and returns the whole thing (so the same remarks applies, in particular: NO ZPVA)

Parameters:

temperature (float) – temperature considered (in Kelvin)

Returns:

(H_t, H_r, H_v) (as a tuple, in Hartree)

Return type:

tuple

compute_entropy(symmetry_number, temperature=298.15, pressure=101325.0)

Compute the value of the different contribution (translation, rotation, vibration) to the entropy, IN HARTREE / KELVIN !

Note

  • Assume that the rotational partition function is equal to the “high” temperature limit (T>100K?)

  • The vibrational contribution does include ZPVA (?)

Parameters:

  • symmetry_number (int) – symmetry number for the rotation

  • temperature (float) – temperature considered (in Kelvin)

  • pressure (float) – pressure (in pascal)

Returns:

(S_e, S_t, S_r, S_v) (as a tuple, in Hartree / Kelvin)

Return type:

tuple

compute_gibbs_free_energy(symmetry_number, temperature=298.15, pressure=101325.0)

Compute the value of the different contribution (translation, rotation, vibration) to gibbs free energy, IN HARTREE !

Note

Since \(G=H-TS\), both enthalpy and entropy are computed, and that is it. No ZPVA, no electronic energy (as usual)

Parameters:

  • symmetry_number (int) – symmetry number for the rotation

  • temperature (float) – temperature considered (in Kelvin)

  • pressure (float) – pressure (in pascal)

Returns:

(G_e, G_t, G_r, G_v) (as a tuple, in Hartree)

Return type:

tuple

compute_internal_energy(temperature=298.15)

Compute the value of the different contribution (translation, rotation, vibration) to the internal energy, IN HARTREE !

Note

  • No electronic contribution

  • Assume that the rotational partition function is equal to the “high” temperature limit (T>100K?)

  • The vibrational contribution does not include ZPVA

Parameters:

temperature (float) – temperature considered (in Kelvin)

Returns:

(U_t, U_r, U_v) (as a tuple, in Hartree)

Return type:

tuple

compute_partition_functions(symmetry_number, temperature=298.15, pressure=101325.0)

Compute the value of the different partition functions

Note

  • Assume that the electronic partition is equal to the multiplicity (T<1000K?)

  • Assume that the rotational partition function is equal to the “high” temperature limit (T>100K?)

  • Assume that the first vibrational level is the zero energy (the “V=0” version of Gaussian)

Parameters:

  • symmetry_number (int) – symmetry number for the rotation

  • temperature (float) – temperature considered (in Kelvin)

  • pressure (float) – pressure (in pascal)

Returns:

(omega_e, omega_t, omega_r, omega_v) (as a tuple)

Return type:

tuple

compute_zpva()

Compute the ZPVA (zero point vibrational averaging) energy from frequencies, but WITHOUT the electronic (SCF) energy

Return type:

float

from_cartesian_hessian(cartesian_hessian)

From a given carthesian hessian, create mass-weigted ones, and extracts frequencies and displacements.

Parameters:

cartesian_hessian (numpy.ndarray) – the cartesian hessian

output_displacements()

Gaussian-like formatting of displacements. Frequencies in \(cm^{-1}\), reduced masses in AMU, displacement in \(\text{Å}\,AMU^{-1/2}\).

Return type:

str

static rotational_temperature(inertia_moment)

Convert an inertia moment (in AMU*angstrom**2) to a rotational temperature (K)

Parameters:

inertia_moment (float) – inertia moment

Return type:

float

static vibrational_temperature(vibrational_energy, scale=1.0)

Convert a vibration energy (hartree) to a vibrational temperature (K)

Parameters:

  • vibrational_energy (float) – the vibrational energy (in Hartree)

  • scale (float) – scale of the vibrational energyµ

Return type:

float

qcip_tools.derivatives_g.REPRESENTATIONS = {'cubic_FF': 'GGG', 'gradient': 'G', 'hessian': 'GG', 'projected gradient': 'N', 'projected hessian': 'NN'}

Correspondence between a name and a representation

qcip_tools.derivatives_g.SIMPLIFIED_NAMES = {'cubic_FF': 'Fijk', 'energy': 'E', 'gradient': 'Fi', 'hessian': 'Hij', 'projected gradient': 'Fi(N)', 'projected hessian': 'Hij(N)'}

List of all simplified names