Tools for molecular vibrations analysis
Project description
Theovib
:construction: Version 0.0.3 :construction:
A package for molecular vibrations analysis. This initial version contains functions to solve the vibrational problem starting from the Hessian matrix or from a 3D-Hessian construct from Interacting Quantum Atoms energy decomposition scheme. Infrared intensities are obtained from atomic charges and dipoles obtained by AIMAll.
Contents
Installation
To install, just use pip:
$ pip install theovib
Usage
You can perform a vibrational analysis for water using the files in the "example" folder.
Vibrational analysis of the Water molecule
In ./example/h2o you will find all single point calculations needed to generate the 3D-Hessian matrix.
The Cartesian coordinates of the water molecule are:
Atom(Label) | X | Y | Z |
---|---|---|---|
O(1) | 0.000000 | 0.000000 | 0.004316 |
H(2) | 0.000000 | -0.763369 | -0.580667 |
H(3) | 0.000000 | 0.763369 | -0.580667 |
The internal coordinates are defined as:
- bond between atoms 1 and 2;
- bond between atoms 1 and 3;
- angle defined by atoms 1, 2 and 3.
We start writing the input file:
MOLECULE:
h2o
FOLDER:
h2o
DELTA:
0.05
BOND:
1 2
1 3
---
ANGLE:
2 1 3
---
DELTA is the displacement that generated the non-equilibrium geometries.
1. Reading the input file:
Import modules:
from theovib.molecule import *
from theovib.internal import *
from theovib.matrices import *
from theovib.ir import *
from theovib.input import *
Use the Input class:
input_data = Input.read_text('input.txt')
Use the Molecule class and its methods to get and store data:
molecule = Molecule.read_gaussian(input_data.folder + '/EQ.com')
molecule.energy = get_energy_from_wfn(input_data.folder +'/EQ.wfn')
molecule.iqa_energy = get_IQA(input_data.folder +'/EQ_atomicfiles', molecule.atoms)
2. Construct the B matrix:
Initialize and construct the B matrix:
b_matrix = []
for coord in input_data.bond:
b_matrix.append(bond(molecule.positions, coord[0], coord[1]))
for coord in input_data.angle:
b_matrix.append(angle(molecule.positions, coord[0], coord[1], coord[2]))
molecule.b_matrix = np.array(b_matrix)
3. Obtain the Hessian and 3D Hessian matrix:
molecule.hessian, molecule.iqa_hessian, errors = hessian_from_iqa(molecule.atoms, input_data.delta, input_data.folder)
The Hessian and 3D Hessian are numerrically generated using IQA contributions from AIMAll outputs
3. Calculate the normal coordinates and infrared intensities:
molecule.normal_coordinates, molecule.freq, molecule.iqa_freq, molecule.iqa_terms = normal_modes(molecule.atoms, molecule.iqa_hessian)
molecule.int, molecule.c_tensors, molecule.ct_tensors, molecule.dp_tensors = intensities(molecule.atoms, molecule.positions, molecule.normal_coordinates, input_data.folder, input_data.delta)
4. Convert to internal coordinates to obtain the force constants:
molecule.internal_hessian, molecule.iqa_forces = convert_to_internal(molecule.atoms, molecule.b_matrix, molecule.iqa_hessian)
Theory
1. Force constant in internal coordinates:
In order to calculate, and decompose the force constants into its IQA components, one needs first to convert the Cartesian Hessian into the the Wilson's $\mathbf{F}$ matrix, that contains the force constants, and their interactions, expressed in internal coordinates. To do this, it is necessary to define the $\mathbf{B}$ matrix, that converts the $N\times 3N$ Cartesian coordinates matrix, $\mathbf{X}$ in the internal coordinates matrix $\mathbf{R}$:
$$ \mathbf{B}\mathbf{X} = \mathbf{R} $$
The process of setting-up the $\mathbf{B}$ matrix is tedious and can be found in the literature. We start by calculating the pseudo-inverse of $\mathbf{B}$ by:
$$ \mathbf{B}^{-1} =\mathbf{M}^{2} \mathbf{B}^{\dagger}\mathbf{G}^{-1} $$
$\mathbf{G}$ contains the inverse of the kinetic energy terms and its inverse, $\mathbf{G}^{-1}$, is given by:
$$\mathbf{G}^{-1} =\mathbf{D} \mathbf{\Phi}^{-1}\mathbf{D}^{\dagger}$$
where $\mathbf{D}$ and $\mathbf{\Phi}$ are, respectively, the eigenvectors and the diagonal eigenvalues matrices of $\mathbf{G}$. The force constants in internal coordinates are then obtained as follows:
$$\mathbf{F} ={\mathbf{B}^{\dagger}}^{-1} \mathbf{H}\mathbf{B}^{-1}$$
The decomposition of the force constants into the IQA contributions is done using:
$$\left[\sum_{k=1}^{N^2} \mathbf{F^{IQA}k}\right] ={\mathbf{B}^{\dagger}}^{-1} \left[\sum{k=1}^{N^2} \mathbf{H^{IQA}_k}\right] \mathbf{B}^{-1}$$
Each $F^{IQA}_k$ is a matrix containing the contribution of the $k^{th}$ IQA term of $\mathbf{F}$.
2. Infrared intensities:
The infrared intensity, $A_k$, of normal mode $Q_k$ is defined as:
$$A_k = \frac{N_A \pi}{3c^2}\left(\frac{\mathrm{d}\vec{p}}{\mathrm{d}Q_k}\right)^2 = \frac{N_A \pi}{3c^2} \sum_{i=1}^3 \left( \frac{\partial\vec{p}}{\partial \sigma_i} \cdot \frac{\partial\sigma_i}{\partial Q_k}\right)^2$$
where $N_A$ is the Avogadro's constant, $c$ is the speed of light and $\vec{p}$ is the molecular dipole moment. The derivative $\frac{\partial\sigma_i}{\partial Q_k}$ is a element of the $k^{th}$ column of $\mathbf{L}$.
$$\mathbf{L} = \mathbf{M}\mathbf{A}$$
Let $\mathbf{T}$ be a block-diagonal matrix of dimension $3N \times 3N$ where each $3 \times 3$ block is an atomic polar tensor. The dipole moment derivative of normal mode $Q_k$ is obtained by solving:
$$\frac{\mathrm{d}p}{\mathrm{d}Q_k} = \mathbf{U} \cdot \mathbf{T} \cdot \mathbf{L_k}$$ $\mathbf{L_k}$ is the $k^{th}$ column of $\mathbf{L}$. $\mathbf{U}$ is a $1 \times N$ line vector whose elements equals 1.
Contributing
Interested in contributing? Check out the contributing guidelines. Please note that this project is released with a Code of Conduct. By contributing to this project, you agree to abide by its terms.
Documentation
https://ljduarte.github.io/theovib/
License
theovib
was created by L. J. Duarte. It is licensed under the terms of the MIT license.
Credits
theovib
was created with cookiecutter
and the py-pkgs-cookiecutter
template.
Project details
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.