Commit aaec79ae authored by psd's avatar psd
Browse files

Exercise 8

parent b95a4bfe
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from ase import Atoms
from ase.io import read,write
from ase.visualize import view
import nglview
```
%% Cell type:markdown id: tags:
## Tasks
### Task 1
Let's begin with some theory. The energy of system of atoms is:
$H = \sum_i p_i/2m_i + VR_i$.
Expanding $V$ around equilibrium ($\sum_i\partial{V}/\partial{R_i}=0$):
$H = V_0 + \frac{1}{2}\sum_{ij} \frac{\partial{V}}{\partial{R_i}\partial{R_j}}\Delta R_i\Delta R_j$
which, in matrix form reads:
$H = V_0 + \frac{1}{2}\{P_{3N}\}\begin{bmatrix}{M_i}^{-1} &&\\ && \\ &&\\ \end{bmatrix}\{P_{3N}\}^T +
\frac{1}{2}\{\Delta R_{3N}\}\begin{bmatrix}{k_{ij}} &&\\ && \\ &&\\ \end{bmatrix}\{\Delta R_{3N}\}^T$
$H = V_0 + \frac{1}{2}\{P_i(\sqrt{M_i})^{-1}\}\{P_i(\sqrt{M_i})^{-1}\}^T +
\frac{1}{2}\{\Delta R_i(\sqrt{M_i})^{-1}\}\underbrace{\begin{bmatrix}{k_{ij}(\sqrt{M_iM_j})^{-1}} &&\\ && \\ &&\\ \end{bmatrix}}_{\text{Hessian matrix}}\{\Delta R_i(\sqrt{M_i})^{-1}\}^T$
In the basis basis states of the Hessian matrix, the above equation can be written:
$H = V_0 + \frac{1}{2}\{\pi_i\}\{\pi_i\}^T +
\frac{1}{2}\{\delta\rho_i\}\begin{bmatrix}{\omega_i} &&\\ &\ddots& \\ &&\\ \end{bmatrix}\{\delta\rho_i\}^T$
where ${\delta\rho_i},{\pi_i}$ are the mass normalized position and momentum vectors. The above equation corresponds to a set of **uncoupled harmonic oscillators**:
$\boxed{H = V_0 + \frac{1}{2}\sum_{\alpha} \Big( \frac{\pi_{\alpha}^2}{2} + \frac{\omega_{\alpha}2}{2}\rho_{\alpha}^2 \Big)}$
The **molden** file contains the eigenvalues and eigenvectors of the Hessian matrix. The eigenvectros corresponds to the atomic displacements and the eigenvalues are the vibrational frequencies. Your task is to devise a method to display the vibrational modes, making use of the *vibr_displacements* and *atoms*(ase.Atoms object) returned by the function *read_molden*.
For the purpose, complete the skeleton function *get_trajectory*, which is initialized in a cell below, and make use of the ase.Atoms class, which is initialized with a list of atomic elements and their respective coordinates.
**Hint**: Write the equation of motion (Lagrangian) for the equation above and find the time evolution of the atoms for a given normal mode.
%% Cell type:markdown id: tags:
### Task 2
- Compare the vibrational frequencies of methanol with experiments (see paper) and the one of benzene with literature on the internet.
- Which kind of modes will correspond to stretching of CH and CC bonds?
- Try to animate some frequencies, and report the kind of mode corresponding to all peaks.
- In the methanol case, you can compare the result you obtained with the one with better basis set and convergence.
- Examine the differences between the file vib.c6h6.inp and the vib.c6h6.ref, and the difference in spectra. Discuss.
%% Cell type:markdown id: tags:
## Funcitons:
**read_molden(file):**
input:
file - molden filename
return:
atoms - ase.Atom object
frequency - vibrational frequencies
vibr_displacements - vibrational atomic displacements in Angstrom
**get_trajectory(mode):**
input:
mode - mode number
return:
trajectory - trajectory of atomic displacements for specified mode
%% Cell type:code id: tags:
``` python
def read_molden(file):
with open(file) as f:
data = f.readlines()
freq = []
info_atoms = [] # element, x, y, z
vibr_displacements = [] # [vibration_nr][coord]
inten = []
section = ''
b2A=0.52917721067
# Parse the datafile
for line in data:
line = line.strip()
# Are we entering into a new section?
if line[0] == '[':
section = line.strip('[]').lower()
continue
if section == 'freq':
freq.append(float(line))
if section == 'fr-coord':
el, x, y, z = line.split()
info_atoms.append([el, float(x)*b2A, float(y)*b2A, float(z)*b2A])
if section == 'fr-norm-coord':
if line.startswith('vibration'):
vibr_displacements.append([])
continue
coords = [float(x) for x in line.split()]
vibr_displacements[-1].append(coords)
if section == 'int':
inten.append(float(line))
vibr_displacements = np.asanyarray(vibr_displacements)
info_atoms = np.asanyarray(info_atoms)
atoms = Atoms(info_atoms[:,0], info_atoms[:,1:4])
return atoms, freq, vibr_displacements
```
%% Cell type:code id: tags:
``` python
def get_trajectory(mode):
'''
TODO: Solve the equation of motion and complete the missing ___
'''
trajectory = []
for time in time_arr:
vibr_atoms = Atoms(a.get_chemical_symbols(), a.positions+'''___''')
trajectory.append(vibr_atoms)
return trajectory
```
%% Cell type:markdown id: tags:
## Read molden file
%% Cell type:code id: tags:
``` python
file = "./C6H6-VIBRATIONS-1.ref.mol"
a, freq, vibr_displacements = read_molden(file)
```
%% Cell type:markdown id: tags:
## View atoms at equilibrium
%% Cell type:code id: tags:
``` python
nglview.show_ase(a)
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
## Plot spectrum
%% Cell type:code id: tags:
``` python
plt.figure(figsize=(10,5))
plt.vlines(freq,np.zeros(len(freq)),np.ones(len(freq)))
plt.hlines(0,*plt.xlim())
plt.xlabel('Energy [cm$^{-1}$]',fontsize=12)
```
%%%% Output: execute_result
Text(0.5, 0, 'Energy [cm$^{-1}$]')
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAlYAAAFKCAYAAADfWRFiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFhlJREFUeJzt3X+wZnddH/D3h10CDiALZFXIbthQQ8cMpZCukRkc3ArCJtZEW7DJ1ApKzdQaqaMyjaVNadrOVBy1g02lQShIlRB+ujDLRBQYaEeSbMoSCWnKJhCyJpLlVzAghNBP/3jOpg+Xu3uf3Xzv3ns3r9fMM/ec7/ne83yeT848+84553ludXcAAHjwHrbWBQAAnCwEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBNq/VE5966qm9Y8eOtXp6AICF3XDDDZ/r7q0rzVuzYLVjx47s27dvrZ4eAGBhVXX7IvNcCgQAGESwAgAYRLACABhEsAIAGESwAgAYRLACABhEsAIAGGTFYFVVr6+qu6vq40fYXlX16qo6UFU3VtXZ48sEAFj/Fjlj9YYku4+y/dwkZ06Pi5P87oMvCwBg41kxWHX3h5J84ShTLkjy+z3zkSRbquqJowoEANgoRtxjdVqSO+bWD05jAAAPKSOCVS0z1stOrLq4qvZV1b5Dhw4NeOq1sWvXrmzZsiW7du3aMHXs2rXriPOOtu3B2LJlS7Zs2bJq+ztc99H6MPq/1XyvRr++Y6lh9PG3Wsf0WvUIWBuH30sOv5/ML6/1v5knyohgdTDJ9rn1bUnuXG5id1/Z3Tu7e+fWrSv+gWgAgA1lRLDak+Snp08HPivJPd1914D9AgBsKJtXmlBVb06yK8mpVXUwyb9J8vAk6e7XJNmb5LwkB5J8NcnPrFaxAADr2YrBqrsvWmF7J/mFYRUBAGxQvnkdAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGCQhYJVVe2uqluq6kBVXbrM9tOr6gNV9dGqurGqzhtfKgDA+rZisKqqTUmuSHJukrOSXFRVZy2Z9q+SXN3dz0xyYZL/MrpQAID1bpEzVuckOdDdt3X3fUmuSnLBkjmd5Dun5ccmuXNciQAAG8PmBeacluSOufWDSX5gyZxXJvnjqvrFJI9K8rwh1QEAbCCLnLGqZcZ6yfpFSd7Q3duSnJfkTVX1bfuuqoural9V7Tt06NCxVwsAsI4tEqwOJtk+t74t336p76VJrk6S7v6zJI9McurSHXX3ld29s7t3bt269fgqBgBYpxYJVtcnObOqzqiqUzK7OX3PkjmfSfLcJKmq78ssWDklBQA8pKwYrLr7/iSXJLkmyc2Zffrvpqq6vKrOn6b9SpKfq6qPJXlzkpd099LLhQAAJ7VFbl5Pd+9NsnfJ2GVzy59I8uyxpQEAbCy+eR0AYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgEMEKAGAQwQoAYBDBCgBgkIWCVVXtrqpbqupAVV16hDk/WVWfqKqbquoPx5YJALD+bV5pQlVtSnJFkh9JcjDJ9VW1p7s/MTfnzCS/luTZ3f3Fqvqu1SoYAGC9WuSM1TlJDnT3bd19X5KrklywZM7PJbmiu7+YJN1999gyAQDWv0WC1WlJ7phbPziNzXtqkqdW1f+sqo9U1e5RBQIAbBQrXgpMUsuM9TL7OTPJriTbkny4qp7W3V/6lh1VXZzk4iQ5/fTTj7lYAID1bJEzVgeTbJ9b35bkzmXm/FF3f6O7P5XklsyC1rfo7iu7e2d379y6devx1gwAsC4tEqyuT3JmVZ1RVackuTDJniVz3pXk7yZJVZ2a2aXB20YWCgCw3q0YrLr7/iSXJLkmyc1Jru7um6rq8qo6f5p2TZLPV9Unknwgycu7+/OrVTQAwHq0yD1W6e69SfYuGbtsbrmT/PL0AAB4SPLN6wAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgAAgywUrKpqd1XdUlUHqurSo8x7YVV1Ve0cVyIAwMawYrCqqk1JrkhybpKzklxUVWctM+8xSV6W5NrRRQIAbASLnLE6J8mB7r6tu+9LclWSC5aZ9++SvCrJ1wbWBwCwYSwSrE5Lcsfc+sFp7AFV9cwk27v7PQNrAwDYUBYJVrXMWD+wsephSX47ya+suKOqi6tqX1XtO3To0OJVAgBsAIsEq4NJts+tb0ty59z6Y5I8LckHq+rTSZ6VZM9yN7B395XdvbO7d27duvX4qwYAWIcWCVbXJzmzqs6oqlOSXJhkz+GN3X1Pd5/a3Tu6e0eSjyQ5v7v3rUrFAADr1IrBqrvvT3JJkmuS3Jzk6u6+qaour6rzV7tAAICNYvMik7p7b5K9S8YuO8LcXQ++LACAjcc3rwMADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMIlgBAAwiWAEADCJYAQAMslCwqqrdVXVLVR2oqkuX2f7LVfWJqrqxqv60qp48vlQAgPVtxWBVVZuSXJHk3CRnJbmoqs5aMu2jSXZ299OTvC3Jq0YXCgCw3i1yxuqcJAe6+7buvi/JVUkumJ/Q3R/o7q9Oqx9Jsm1smQAA698iweq0JHfMrR+cxo7kpUneu9yGqrq4qvZV1b5Dhw4tXiUAwAawSLCqZcZ62YlVP5VkZ5LfWG57d1/Z3Tu7e+fWrVsXrxIAYAPYvMCcg0m2z61vS3Ln0klV9bwkr0jyQ9399THlAQBsHIucsbo+yZlVdUZVnZLkwiR75idU1TOT/Nck53f33ePLBABY/1YMVt19f5JLklyT5OYkV3f3TVV1eVWdP037jSSPTvLWqtpfVXuOsDsAgJPWIpcC0917k+xdMnbZ3PLzBtcFALDh+OZ1AIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBFgpWVbW7qm6pqgNVdeky2x9RVW+Ztl9bVTtGFwoAsN6tGKyqalOSK5Kcm+SsJBdV1VlLpr00yRe7+3uT/HaSXx9dKADAerfIGatzkhzo7tu6+74kVyW5YMmcC5K8cVp+W5LnVlWNKxMAYP3bvMCc05LcMbd+MMkPHGlOd99fVfckeUKSz40o8njt2rVrVfa7f//+3Hvvvdm/f/+qPcfoOvbv359k+Z4cbduDce+99w7d79L9Ha778Lbl+jD6v9V8r0a/vmOpYfTxt1rH9Fr1CFgbh99LlltOVu+94IMf/OCq7Pd4VHcffULVi5K8oLv/ybT+j5Oc092/ODfnpmnOwWn91mnO55fs6+IkFyfJ6aef/nduv/32ka/l23gzB4CT34kIVlV1Q3fvXGneImesDibZPre+LcmdR5hzsKo2J3lski8s3VF3X5nkyiTZuXPn0RPdAOspwQIAJ79F7rG6PsmZVXVGVZ2S5MIke5bM2ZPkxdPyC5O8v1c6FQYAcJJZ8YzVdM/UJUmuSbIpyeu7+6aqujzJvu7ek+R1Sd5UVQcyO1N14WoWDQCwHi1yKTDdvTfJ3iVjl80tfy3Ji8aWBgCwsfjmdQCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQVb8W4Gr9sRVh5Ks7h8LXP9OzRr/oeqTkJ6Op6fj6elY+jmenn67J3f31pUmrVmwIqmqfYv8QUcWp6fj6el4ejqWfo6np8fPpUAAgEEEKwCAQQSrtXXlWhdwEtLT8fR0PD0dSz/H09Pj5B4rAIBBnLECABhEsFplVfXpqvrzqtpfVfumscdX1fuq6pPTz8dN41VVr66qA1V1Y1WdvbbVrw9V9fqquruqPj43dsw9rKoXT/M/WVUvXovXsh4coZ+vrKq/mI7T/VV13ty2X5v6eUtVvWBufPc0dqCqLj3Rr2M9qartVfWBqrq5qm6qqn8+jTtOj8NR+uk4PU5V9ciquq6qPjb19N9O42dU1bXT8faWqjplGn/EtH5g2r5jbl/L9ppJd3us4iPJp5OcumTsVUkunZYvTfLr0/J5Sd6bpJI8K8m1a13/engkeU6Ss5N8/Hh7mOTxSW6bfj5uWn7cWr+2ddTPVyb51WXmnpXkY0kekeSMJLcm2TQ9bk3ylCSnTHPOWuvXtoY9fWKSs6flxyT5P1PvHKdj++k4Pf6eVpJHT8sPT3LtdOxdneTCafw1SX5+Wv5nSV4zLV+Y5C1H6/Vav7719HDGam1ckOSN0/Ibk/z43Pjv98xHkmypqieuRYHrSXd/KMkXlgwfaw9fkOR93f2F7v5ikvcl2b361a8/R+jnkVyQ5Kru/np3fyrJgSTnTI8D3X1bd9+X5Kpp7kNSd9/V3f9rWv6rJDcnOS2O0+NylH4eieN0BdOxdu+0+vDp0Ul+OMnbpvGlx+jhY/dtSZ5bVZUj95qJYLX6OskfV9UNVXXxNPbd3X1XMnsDSfJd0/hpSe6Y+92DOfqbyUPZsfZQb1d2yXRZ6vWHL1lFP4/ZdMnkmZmdEXCcPkhL+pk4To9bVW2qqv1J7s4stN+a5Evdff80Zb4/D/Ru2n5PkidET1ckWK2+Z3f32UnOTfILVfWco8ytZcZ8bPPYHKmHent0v5vkbyR5RpK7kvzmNK6fx6CqHp3k7Ul+qbu/fLSpy4zp6xLL9NNx+iB09ze7+xlJtmV2lun7lps2/dTT4yRYrbLuvnP6eXeSd2Z2MH/28CW+6efd0/SDSbbP/fq2JHeeuGo3lGPtod4eRXd/dnrT/b9JXpv/f2pfPxdUVQ/PLAT8QXe/Yxp2nB6n5frpOB2ju7+U5IOZ3WO1pao2T5vm+/NA76btj83sFgI9XYFgtYqq6lFV9ZjDy0men+TjSfYkOfxpnxcn+aNpeU+Sn54+MfSsJPccvozAtznWHl6T5PlV9bjp8sHzpzHywD/6h/1EZsdpMuvnhdMnhM5IcmaS65Jcn+TM6RNFp2R2c+ueE1nzejLde/K6JDd392/NbXKcHocj9dNxevyqamtVbZmWvyPJ8zK7d+0DSV44TVt6jB4+dl+Y5P3d3Tlyrzlsre+eP5kfmX0S5WPT46Ykr5jGn5DkT5N8cvr5+Gm8klyR2XXvP0+yc61fw3p4JHlzZqf9v5HZ/y299Hh6mORnM7vR8kCSn1nr17XO+vmmqV83ZvbG+cS5+a+Y+nlLknPnxs/L7NNatx4+th+qjyQ/mNnlkBuT7J8e5zlOh/fTcXr8PX16ko9Ovft4ksum8adkFowOJHlrkkdM44+c1g9M25+yUq89Zg/fvA4AMIhLgQAAgwhWAACDCFYAAIMIVgAAgwhWAACDCFYAAIMIVgCDVdVjq+q6qrq3qp621vUAJ45gBTDeV5P8aJK3rXUhwIklWAEntarqqvpKVf2HE/Wc3f2N7j60TC3vr6qvVdX/OFG1ACeWYAUspKo+XVV/PV3eOvz4z2td14L+dne/Yq2L6O4fTvJP17oOYPVsXnkKwAN+rLv/ZLV2XlWbu/v+1dr/SFX1PVn+Ut8Lu/svT3Q9wPrgjBXwoE1ns361qm6sqnuq6i1V9ci57U+qqrdX1aGq+lRVvWzJ7/6LqroxyVeqanNVnV1VH62qv6qqt077+/fT/JdX1duXPP/vVNV/WrDW7VX1jqmWz8+fdZtqefn0Or5SVa+rqu+uqvdOtfxJVT0uSbr7L7v7B5d5CFXwECZYAaP8ZJLdSc5I8vQkL0mSqnpYkncn+ViS05I8N8kvVdUL5n73osxu9t6S2fvSO5O8Icnjk7w5yU/Mzf3vSXZX1ZZp/5uT/MMkb1qpwKralOQ9SW5PsmOq56ol0/5Bkh9J8tQkP5bkvUn+ZZJTp9pelgVU1d4kz0/y2qp6ySK/A2x8LgUCx+JdVTV/qe7l3f3aafnV3X1nklTVu5M8Yxr//iRbu/vyaf22qnptkguTXDP3u3dMv/uczN6bXt3dneQdVXXd4Sfs7ruq6kNJXpTktZmFuc919w0L1H9OkidNdR9+HUtvJP+d7v7sVMuHk9zd3R+d1t+ZWTBcUXeft8g84OQiWAHH4sePco/V/CWwr2YWYJLkyUmeVFVfmtu+KcmH59bvmFt+UpK/mELVctuT5I1Jfj6zYPVTWeBs1WR7kttXuI/rs3PLf73M+qMXfC7gIcilQGC13ZHkU929Ze7xmCVndOZD1F1JTquqmhvbvmSf70ry9OnLN/9ekj84hlpOny4fAgwnWAGr7bokX55uUP+OqtpUVU+rqu8/wvw/S/LNJJdMN7JfkNklvAd099cy+0TeHya5rrs/cwy13JXkP1bVo6rqkVX17ON6VQDLEKyAY/HuJd9j9c6VfqG7v5nZTeDPSPKpJJ9L8ntJHnuE+fcl+ftJXprkS5ld6ntPkq8vmfrGJH8ri18GnK/le5N8JsnBzG58BxiivvU2BoD1p6quTfKa7v5vc2OnJ/nfSb6nu798lN/9Wmah7NXd/a9XvdijqKr3JXlWZmfZFroJHthYBCtg3amqH0pyS2Znt/5RktckeUp33zVtf1iS30rynd39s2tWKMASbuAE1qO/meTqzD6Bd2tm32Z+OFQ9KrNP6t2e2VctAKwbzlgBAAzi5nUAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEEEKwCAQQQrAIBBBCsAgEH+HzM1wXDKIRITAAAAAElFTkSuQmCC)
%% Cell type:markdown id: tags:
## Visualize vibrational mode
%% Cell type:code id: tags:
``` python
mode = 3
trajectory = get_trajectory(mode)
nglv = nglview.show_asetraj(trajectory, gui=True)
nglv
```
%%%% Output: display_data
%%%% Output: display_data
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from ase import Atoms
from ase.io import read,write
from ase.visualize import view
import nglview
```
%% Cell type:markdown id: tags:
## Tasks
### Task 1
Let's begin with some theory. The energy of system of atoms is:
$H = \sum_i p_i/2m_i + VR_i$.
Expanding $V$ around equilibrium ($\sum_i\partial{V}/\partial{R_i}=0$):
$H = V_0 + \frac{1}{2}\sum_{ij} \frac{\partial{V}}{\partial{R_i}\partial{R_j}}\Delta R_i\Delta R_j$
which, in matrix form reads:
$H = V_0 + \frac{1}{2}\{P_{3N}\}\begin{bmatrix}{M_i}^{-1} &&\\ && \\ &&\\ \end{bmatrix}\{P_{3N}\}^T +
\frac{1}{2}\{\Delta R_{3N}\}\begin{bmatrix}{k_{ij}} &&\\ && \\ &&\\ \end{bmatrix}\{\Delta R_{3N}\}^T$
$H = V_0 + \frac{1}{2}\{P_i(\sqrt{M_i})^{-1}\}\{P_i(\sqrt{M_i})^{-1}\}^T +
\frac{1}{2}\{\Delta R_i(\sqrt{M_i})^{-1}\}\underbrace{\begin{bmatrix}{k_{ij}(\sqrt{M_iM_j})^{-1}} &&\\ && \\ &&\\ \end{bmatrix}}_{\text{Hessian matrix}}\{\Delta R_i(\sqrt{M_i})^{-1}\}^T$
In the basis basis states of the Hessian matrix, the above equation can be written:
$H = V_0 + \frac{1}{2}\{\pi_i\}\{\pi_i\}^T +
\frac{1}{2}\{\delta\rho_i\}\begin{bmatrix}{\omega_i} &&\\ &\ddots& \\ &&\\ \end{bmatrix}\{\delta\rho_i\}^T$
where ${\delta\rho_i},{\pi_i}$ are the mass normalized position and momentum vectors. The above equation corresponds to a set of **coupled harmonic hoscillators**:
$\boxed{H = V_0 + \frac{1}{2}\sum_{\alpha} \Big( \frac{\pi_{\alpha}^2}{2} + \frac{\omega_{\alpha}2}{2}\rho_{\alpha}^2 \Big)}$
The **molden** file contains the eigenvalues and eigenvectors of the Hessian matrix. The eigenvectros corresponds to the atomic displacements and the eigenvalues are the vibrational frequencies. Your task is to devise a method to display the vibrational modes, making use of the *vibr_displacements* and *atoms*(ase.Atoms object) return by the function *read_molden*.
For the purpose, complete the skeleton function *get_trajectory*, which is initialized in a cell below, and make use of the ase.Atoms class, which is initialized with a list of atomic elements and their respective coordinates.
**Hint**: Write the equation of motion (Lagrangian) for the equation above and find the time evolution of the atoms for a given normal mode.
%% Cell type:markdown id: tags:
### Task 2
- Compare the vibrational frequencies of methanol with experiments (see paper) and the one of benzene with literature on the internet.
- Which kind of modes will correspond to stretching of CH and CC bonds?
- Try to animate some frequencies, and report the kind of mode corresponding to all peaks.
- In the methanol case, you can compare the result you obtained with the one with better basis set and convergence.
- Examine the differences between the file vib.c6h6.inp and the vib.c6h6.ref, and the difference in spectra. Discuss.
%% Cell type:markdown id: tags:
## Funcitons:
**read_molden(file):**
input:
file - molden filename
return:
atoms - ase.Atom object
frequency - vibrational frequencies
vibr_displacements - vibrational atomic displacements in Angstrom
**get_trajectory(mode):**
input:
mode - mode number
return:
trajectory - trajectory of atomic displacements for specified mode
%% Cell type:code id: tags:
``` python
def read_molden(file):
with open(file) as f:
data = f.readlines()
freq = []
info_atoms = [] # element, x, y, z
vibr_displacements = [] # [vibration_nr][coord]
inten = []
section = ''
b2A=0.52917721067
# Parse the datafile
for line in data:
line = line.strip()
# Are we entering into a new section?
if line[0] == '[':
section = line.strip('[]').lower()
continue
if section == 'freq':
freq.append(float(line))
if section == 'fr-coord':
el, x, y, z = line.split()
info_atoms.append([el, float(x)*b2A, float(y)*b2A, float(z)*b2A])
if section == 'fr-norm-coord':
if line.startswith('vibration'):
vibr_displacements.append([])
continue
coords = [float(x) for x in line.split()]
vibr_displacements[-1].append(coords)
if section == 'int':
inten.append(float(line))
vibr_displacements = np.asanyarray(vibr_displacements)
embed()
info_atoms = np.asanyarray(info_atoms)
atoms = Atoms(info_atoms[:,0], info_atoms[:,1:4])
return atoms, freq, vibr_displacements
```
%% Cell type:code id: tags:
``` python
def get_trajectory(mode):
enhance_disp = 2.0
time_arr = np.linspace(0.0, 2*np.pi, 20)
trajectory = []
for time in time_arr:
#eq_atoms = ase.Atoms(len(chem_symbols)*['O'], eq_coords)
vibr_atoms = Atoms(a.get_chemical_symbols(), a.positions+enhance_disp*np.sin(time)*vibr_displacements[mode])
trajectory.append(vibr_atoms)
return trajectory
```
%% Cell type:markdown id: tags:
## Read molden file
%% Cell type:code id: tags:
``` python
file = "./C6H6-VIBRATIONS-1.ref.mol"
a, vibr_displacements = read_molden(file)
```
%%%% Output: stream
Python 3.7.1 (default, Dec 14 2018, 19:28:38)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.2.0 -- An enhanced Interactive Python. Type '?' for help.
In [1]: info_atoms
Out[1]:
[['C', 2.222544284814e-05, 1.396857970284175, 0.0],
['H', -1.2171075845410001e-05, 2.487163053250008, -0.0],
['C', -1.210924148834321, 0.6987107720067693, -0.0],
['H', -2.1548170103291895, 1.2443078223466488, 0.0],
['C', -1.2109310281380596, -0.6986636752350197, 0.0],
['H', -2.154787376405392, -1.2443226393085474, -0.0],
['C', -2.169626563747e-05, -1.396875433132127, -0.0],
['H', 1.2171075845410001e-05, -2.4871815744523817, 0.0],
['C', 1.210909861049633, -0.6987023051713985, 0.0],
['H', 2.1548021933672907, -1.244299355511278, 0.0],
['C', 1.2109463742771691, 0.6986721420703903, -0.0],
['H', 2.1548048392533437, 1.2443311061439182, 0.0]]
In [2]: exit
%% Cell type:markdown id: tags:
## View atoms at equilibrium
%% Cell type:code id: tags:
``` python
nglview.show_ase(a)
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
## Plot spectrum
%% Cell type:code id: tags:
``` python
plt.figure(figsize=(10,5))
plt.vlines(freq,np.zeros(len(freq)),np.ones(len(freq)))
plt.hlines(0,*plt.xlim())
plt.xlabel('Energy [cm$^{-1}$]',fontsize=12)
```
%%%% Output: error
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-6-8c677f15f735> in <module>
1 plt.figure(figsize=(10,5))
----> 2 plt.vlines(freq,np.zeros(len(freq)),np.ones(len(freq)))
3 plt.hlines(0,*plt.xlim())
4 plt.xlabel('Energy [cm$^{-1}$]',fontsize=12)
NameError: name 'freq' is not defined
%%%% Output: display_data
%% Cell type:markdown id: tags:
## Visualize vibrational mode
%% Cell type:code id: tags:
``` python
mode = 3
trajectory = get_trajectory(mode)
nglv = nglview.show_asetraj(trajectory, gui=True)
nglv
```
This diff is collapsed.
[Molden Format]
[FREQ]
400.302903
401.555638
596.097032
602.918089
661.189688
706.249913
833.568431
833.756231
949.504594
949.976122
976.963838
985.823710
993.020931
1029.977506
1035.062465
1143.899937
1163.645030
1165.969536
1338.887202
1347.664359
1466.711983
1471.797438
1593.030396
1594.117474
3087.283723
3096.781439
3097.757885
3111.829664
3112.819108
3121.735327
[FR-COORD]
C 0.000042 2.639679 0.000000
H -0.000023 4.700057 -0.000000
C -2.288315 1.320372 -0.000000
H -4.072014 2.351401 0.000000
C -2.288328 -1.320283 0.000000
H -4.071958 -2.351429 -0.000000
C -0.000041 -2.639712 -0.000000
H 0.000023 -4.700092 0.000000
C 2.288288 -1.320356 0.000000
H 4.071986 -2.351385 0.000000
C 2.288357 1.320299 -0.000000
H 4.071991 2.351445 0.000000
[FR-NORM-COORD]
vibration 1
-0.000005 0.000033 0.000502
0.000008 0.000052 0.001115
-0.000006 -0.000005 0.201939
-0.000014 -0.000001 0.456716
0.000029 -0.000003 -0.202439
0.000029 0.000023 -0.457874
-0.000000 -0.000010 0.000502
-0.000007 -0.000001 0.001115
-0.000015 -0.000015 0.201939
-0.000009 0.000004 0.456736
-0.000005 -0.000005 -0.202439
0.000013 -0.000022 -0.457854
vibration 2
-0.000012 0.000028 0.233236
0.000011 0.000041 0.528255
0.000003 -0.000006 -0.117066
-0.000001 -0.000018 -0.264924
0.000024 0.000003 -0.116197
0.000027 0.000024 -0.262997
-0.000004 -0.000007 0.233236
-0.000007 0.000001 0.528231
-0.000005 -0.000015 -0.117065
-0.000000 -0.000012 -0.264935
-0.000010 -0.000006 -0.116197
0.000001 -0.000016 -0.262993
vibration 3
-0.000042 0.357385 -0.000048
0.000111 0.355791 0.000060
0.224264 0.024295 -0.000010
0.058038 -0.261826 -0.000001
0.224428 -0.024164 0.000061
0.058276 0.261838 -0.000101
0.000053 -0.357363 -0.000019
-0.000009 -0.355772 -0.000057
-0.224283 -0.024294 -0.000027
-0.058128 0.261748 0.000037
-0.224421 0.024148 0.000033
-0.058266 -0.261852 0.000170
vibration 4
-0.151123 -0.000115 -0.000032
0.231673 -0.000110 -0.000092
-0.232074 0.220707 0.000016
-0.325667 0.056506 0.000006
0.231918 0.220718 0.000016
0.325563 0.056426 -0.000152
0.151125 0.000103 0.000009
-0.231598 0.000084 -0.000037
0.232078 -0.220697 -0.000015
0.325650 -0.056526 -0.000118
-0.231927 -0.220715 0.000036
-0.325603 -0.056390 0.000042
vibration 5
0.000011 0.000006 -0.034007
0.000015 0.000015 0.406053
-0.000017 0.000015 -0.034214
-0.000033 0.000015 0.407191
-0.000007 -0.000008 -0.034184
-0.000006 -0.000030 0.407242
0.000014 0.000002 -0.034041
0.000007 0.000007 0.406001
0.000009 -0.000003 -0.034180
0.000006 0.000002 0.407241
-0.000009 -0.000011 -0.034216
-0.000007 -0.000026 0.407180
vibration 6
-0.000010 0.000006 0.201113
0.000012 0.000009 0.356838
-0.000011 0.000016 -0.201393
-0.000026 -0.000015 -0.354372
0.000018 0.000010 0.201408
0.000023 -0.000002 0.354173
0.000008 -0.000002 -0.201095
-0.000016 0.000000 -0.357067
0.000008 -0.000018 0.201412
0.000021 0.000012 0.354145
-0.000013 -0.000013 -0.201388
-0.000016 0.000001 -0.354400
vibration 7
-0.000006 0.000021 -0.000410
-0.000022 0.000024 0.002727
0.000002 0.000001 -0.074048
-0.000015 -0.000024 0.495897
0.000005 0.000003 -0.073645
0.000009 0.000001 0.493120
-0.000008 -0.000019 0.000410
-0.000006 -0.000020 -0.002735
0.000004 -0.000014 0.074048
0.000024 0.000005 -0.495875
0.000004 0.000007 0.073647
0.000007 0.000018 -0.493151
vibration 8
-0.000002 0.000013 0.085376
0.000002 0.000012 -0.569602
0.000000 -0.000006 0.042262
-0.000014 -0.000024 -0.284524
0.000029 -0.000002 -0.042966
0.000034 0.000006 0.289248
-0.000001 0.000010 -0.085378
-0.000010 0.000026 0.569633
-0.000020 -0.000008 -0.042260
-0.000023 -0.000004 0.284498
-0.000007 -0.000008 0.042966
0.000008 -0.000026 -0.289260
vibration 9