Vibrational Corrections to Properties

In order to compare experimental calculations with theoretical calculations at a high accuracy, one needs to account for the vibrational motion of the nuclei. With this implementation based on perturbation theory to second order (VPT2), one can account for this vibrational motion and correct the equilibrium property value ad hoc. The vibrational correction to the property \(P\) is given as:

\[\Delta^\text{VPT2}P = -\frac{1}{2}\sum_i\frac{1}{\omega_i}\frac{\partial P}{\partial q_i}\sum_j{k_{ijj}}\left(\nu_j+\frac{1}{2}\right) + \frac{1}{2}\sum_{i}\left(\nu_i+\frac{1}{2}\right)\frac{\partial^2 P}{\partial q_i^2}\]

where \(k_{ijj}\) is the cubic force constant in \(\text{cm}^{-1}\) and \(\omega\) is the harmonic frequency in \(\text{cm}^{-1}\) with all derivatives evaluated at \(q=0\). For a zero-point vibrational correction (ZPVC), \(\nu\) is set to 0. Temperature dependent corrections with a rotational contribution can also be computed using a Boltzmann averaging of states:

\[\Delta^\text{VPT2}P = -\frac{1}{4}\sum_i\frac{1}{\omega_i}\frac{\partial P}{\partial q_i}\left(\sum_j{k_{ijj}}\coth\left(\frac{{h}c\omega_j}{2kT}\right) - \frac{kT}{2\pi c}\frac{1}{\sqrt{hc\omega_i}}\sum_\alpha \frac{a^{\alpha\alpha}_i}{I^e_{\alpha\alpha}}\right)\]
\[+ \frac{1}{4}\sum_{i}\frac{\partial^2 P}{\partial q_i^2}\coth\left(\frac{{h}c\omega_i}{2kT}\right)\]

where \(T\) is the temperature in Kelvin, \(k\) is the Boltzmann constant, \(h\) is the Planck constant, \(I\) is the tensor of the principle axes of inertia and \(a\) is the tensor of linear expansion coefficients of the moment of inertia in the normal coordinates. Cubic force constants, first and second derivatives of properties are derived numerically from displaced geometry calculations using either a 3/5 point stencil or a polynomial fitting. The stepsize used for the numerical analysis can be given as an argument and is unitless (reduced normal coordinates).

The following properties are currently supported (program dependent):

  • NMR shieldings - Gaussian/Dalton

  • Spin-spin coupling constants - Gaussian/Dalton

  • Hyperfine coupling constants - Gaussian

  • Static polarizability/frequency-dependent polarizabilities - Dalton/Gaussian

  • Near static optical rotation/frequency-dependent optical rotations - Dalton/Gaussian

Tips:

  • Performing a geometry optimisation is always advised (unless providing an already optimised geometry).

  • The temperature can be given as an argument for calculations of temperature-dependent corrections.

  • The fitting method (linear/polynomial) can be changed to improve the accuracy of the results.

  • Temperature, fitting method, fitting order/point stencil can be changed, and the script can be rerun interactively without explicitly calling a quantum chemistry program.

  • If the stepsize is changed or the input geometry is changed in any way, then all files need to be deleted.

  • Polynomial fitting results can be inspected by changing the plot_polyfittings variable to “True” in the input file.

  • It is advised to create an output folder for the generated files and run the script from there (as many files can be generated!)

Note, that the quantum chemistry methods used in this tutorial may not be suitable for the particular problem that you want to solve. Moreover, the basis set used here is minimal in order to allow you to progress fast through the tutorial and should under no circumstances be used in production calculations.

Correction to NMR shieldings using a linear fitting

import daltonproject as dp

#  Essential settings
hf = dp.QCMethod(qc_method='HF')
basis = dp.Basis(basis='STO-3G')
molecule = dp.Molecule(input_file='water.xyz')

# Settings if running multiple calculations in parallel using job farming
settings = dp.ComputeSettings(mpi_num_procs=1, jobs_per_node=1)

# optional - optimize the geometry
geo_opt = dp.Property(geometry_optimization=True)
opt_output = dp.dalton.compute(molecule, basis, hf, geo_opt, compute_settings=settings)
molecule.coordinates = opt_output.final_geometry

# Essential settings - compute the reference geometry Hessian
hess = dp.Property(hessian=True)
prop_output = dp.dalton.compute(molecule, basis, hf, hess, compute_settings=settings)
vib_prop = dp.Property(nmr_shieldings=True)

# Essential settings for computing vibrational correction
va_settings = dp.VibAvSettings(molecule=molecule,
                               property_program='dalton',
                               is_mol_linear=False,
                               hessian=prop_output.hessian,
                               property_obj=vib_prop,
                               stepsize=0.05,
                               differentiation_method='linear',
                               temperature=0,
                               linear_point_stencil=3)

# Essential settings - Instances of Molecule class for distorted geometries
molecules = []
for i in va_settings.file_list:
    molecules.append(dp.Molecule(input_file=i))

# Essential settings - Compute Hessians + property at distorted geometries
dalton_results_hess = []
dalton_results_nmr = []
for mol in molecules:
    dalton_results_hess.append(dp.dalton.compute(mol, basis, hf, hess, compute_settings=settings))
    dalton_results_nmr.append(dp.dalton.compute(mol, basis, hf, vib_prop, compute_settings=settings))

# Essential settings - Instance of ComputeVibAvCorrection class containing the correction
result = dp.ComputeVibAvCorrection(dalton_results_hess, dalton_results_nmr, va_settings)

print('\nVibrational Correction to Property:', result.vibrational_corrections)

To start, we will calculate corrections to NMR shieldings for water. Here the cubic force constants and property derivatives will be calculated using a 3-point linear stencil with a stepsize of 0.05. An output file containing all corrected shieldings will be generated.

Correction to spin-spin coupling constants using a polynomial fitting

import daltonproject as dp

# Optional - set isotopes for each atom
atom_isotopes = ['O17', 'H1', 'H1']

#  Essential settings
hf = dp.QCMethod(qc_method='HF')
basis = dp.Basis(basis='STO-3G')
molecule = dp.Molecule(input_file='water.xyz', isotopes=atom_isotopes)

# Compute settings
settings = dp.ComputeSettings(mpi_num_procs=1, jobs_per_node=1)

# optional - optimize the geometry
geo_opt = dp.Property(geometry_optimization=True)
opt_output = dp.dalton.compute(molecule, basis, hf, geo_opt, compute_settings=settings)
molecule.coordinates = opt_output.final_geometry

# Essential settings - compute the reference geometry Hessian
hess = dp.Property(hessian=True)
prop_output = dp.dalton.compute(molecule, basis, hf, hess, compute_settings=settings)
vib_prop = dp.Property(spin_spin_couplings=True)

# Essential settings for computing vibrational correction
va_settings = dp.VibAvSettings(molecule=molecule,
                               property_program='dalton',
                               is_mol_linear=False,
                               hessian=prop_output.hessian,
                               property_obj=vib_prop,
                               stepsize=0.05,
                               differentiation_method='polynomial',
                               temperature=298,
                               polynomial_fitting_order=3,
                               plot_polyfittings=True)

# Essential settings - Instances of Molecule class for distorted geometries
molecules = []
for i in va_settings.file_list:
    molecules.append(dp.Molecule(input_file=i, isotopes=atom_isotopes))

# Essential settings - Compute Hessians + property at distorted geometries
dalton_result_hess = []
dalton_result_nmr = []
for mol in molecules:
    dalton_result_hess.append(dp.dalton.compute(mol, basis, hf, hess, compute_settings=settings))
    dalton_result_nmr.append(dp.dalton.compute(mol, basis, hf, vib_prop, compute_settings=settings))

# Essential settings - Instance of ComputeVibAvCorrection class containing the correction
result = dp.ComputeVibAvCorrection(dalton_result_hess, dalton_result_nmr, va_settings)

print('\nVibrational Correction to Property:', result.vibrational_corrections)

This time we will calculate corrections to the spin-spin coupling constants for water. The cubic force constants and property derivatives will be calculated using a 3rd-order polynomial fitting with a stepsize of 0.05. Plots of the property derivatives fitting can also be generated for inspection to ensure the fitting is resonable. Here we also altered the temperature to 298K to calculate the temperature-dependent corrections. An output file containing all corrected spin-spin coupling constants will be generated.

Graphs of the polynomial derivative fitting for each coupling and associated mode is as shown:

../../_images/fc_corr_water.svg

Correction to frequency-dependent polarizabilities using a linear fitting

import daltonproject as dp

# Optional - set isotope for each atom
# Most abundant isotopes are used by default
atom_isotopes = ['O17', 'H1', 'H1']

#  Essential settings
hf = dp.QCMethod(qc_method='HF', scf_threshold=1e-10)
basis = dp.Basis(basis='STO-3G')
molecule = dp.Molecule(input_file='water.xyz', isotopes=atom_isotopes)

# Compute settings
settings = dp.ComputeSettings(mpi_num_procs=1, jobs_per_node=1)

# optional - optimize the geometry
geo_opt = dp.Property(geometry_optimization=True)
opt_output = dp.dalton.compute(molecule, basis, hf, geo_opt, compute_settings=settings)
molecule.coordinates = opt_output.final_geometry

# Essential settings - compute the reference geometry Hessian
hess = dp.Property(hessian=True)
prop_output = dp.dalton.compute(molecule, basis, hf, hess, compute_settings=settings)
vib_prop = dp.Property(polarizabilities={'frequencies': [0.1, 0.2, 0.3]})

# Essential settings for computing vibrational correction
va_settings = dp.VibAvSettings(
    molecule=molecule,
    property_program='dalton',
    is_mol_linear=False,
    hessian=prop_output.hessian,
    property_obj=vib_prop,
    stepsize=0.05,
    differentiation_method='linear',
    temperature=0,
    # polynomial_fitting_order=3,
    linear_point_stencil=3,
    plot_polyfittings=True)

# Essential settings - Instances of Molecule class for distorted geometries
molecules = []
for i in va_settings.file_list:
    molecules.append(dp.Molecule(input_file=i, isotopes=atom_isotopes))

dalton_result_hess = []
dalton_result_nmr = []
for mol in molecules:
    dalton_result_hess.append(dp.dalton.compute(mol, basis, hf, hess, compute_settings=settings))
    dalton_result_nmr.append(dp.dalton.compute(mol, basis, hf, vib_prop, compute_settings=settings))

# Essential settings - Instance of ComputeVibAvCorrection class containing the correction
result = dp.ComputeVibAvCorrection(dalton_result_hess, dalton_result_nmr, va_settings)

print('\nVibrational Correction to the Property:', result.vibrational_corrections)
# result.vibrational_correction =  [0.126377, 0.132178, 0.15234, 0.199484]

import numpy as np  # noqa

np.testing.assert_allclose(result.vibrational_corrections, [0.126377, 0.132178, 0.15234, 0.199484], atol=1e-6)

This time, we will calculate corrections to the frequency-dependent polarizabilities for water where frequencies are in a.u. Static polarizabilities are also done by default when frequency-dependent polarizabilities are requested.

Optical rotations can be requested by changing

vib_prop = dp.Property(polarizabilities={'frequencies': [0.1, 0.2, 0.3]})

to

vib_prop = dp.Property(optical_rotations={'frequencies': [0.1, 0.2, 0.3]})

Ref: Faber, R.; Kaminsky, J.; Sauer, S. P. A. In Gas Phase NMR, Jackowski, K., Jaszunski, M., Eds.; Royal Society of Chemistry, London: 2016; Chapter 7, pp 218–266