# 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:

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 can also be computed using a Boltzmann averaging of states:

where \(T\) is the temperature in Kelvin, \(k\) is the Boltzmann constant and \(h\) is Planks constant. 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:

NMR shieldings

Spin-spin coupling constants

Polarizabilities

- 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_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)
```

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:

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