"Methods for solving the one dimentional vibrational eigenproblem"
from __future__ import print_function, division, absolute_import
import numpy as np
import pandas as pd
from scipy.constants import value, Boltzmann, Avogadro, Planck, gas_constant
from .io import print_mode_thermo
[docs]def factsqrt(m, n):
"""
Return a factorial like constant
Parameters
----------
m : int
Argument of the series
n : int
Length of the series
Notes
-----
.. math::
f(m, n) = \prod^{n - 1}_{i = 0} \sqrt{m - i}
"""
return np.sqrt(np.prod([m - i for i in range(n)]))
[docs]def get_hamiltonian(rank, freq, mass, coeffs):
"""
Compose the Hamiltonian matrix for the anharmonic oscillator with the
potential described by the sixth order polynomial.
Parameters
----------
rank : int
Rank of the Hamiltonian matrix
freq : float
Fundamental frequency in hartrees
mass : float
Reduced mass of the mode
coeffs : array
A one dimensional array with polynomial coeffients
Notes
-----
.. math::
H_{ij} = \left\langle \Psi_{i} \left| \hat{H} \\right| \Psi_{j} \\right\\rangle
where
.. math::
\hat{H} = -\\frac{\hbar^2}{2}\\frac{\partial^2}{\partial \\boldsymbol{Q}^2}
+ \sum_{\mu=0}^{6}c_{\mu}\\boldsymbol{Q}^{\mu}
and :math:`\Psi_{i}` are the standard harmonic oscillator functions.
"""
Hamil = np.zeros((rank, rank), dtype=float)
# change this to proper value
vk = np.sqrt(1.0 / (mass * freq))
uk = -0.5 * freq
# main diagonal i == j
idx = np.arange(rank)
Hamil[idx, idx] = [
-0.5 * uk * (2 * i + 1.0)
+ coeffs[0]
+ 0.5 * coeffs[2] * vk ** 2 * (2 * i + 1.0)
+ 0.25 * coeffs[4] * vk ** 4 * (6.0 * i ** 2 + 6.0 * i + 3)
+ 0.125
* coeffs[6]
* vk ** 6
* (20.0 * i ** 3 + 30.0 * i ** 2 + 40.0 * i + 15.0)
for i in idx
]
# diagonal wih offset 1, i == j + 1
k = rank - 1
idx = np.arange(k)
Hamil[idx + 1, idx] = [
np.sqrt(2.0 * i)
* (
0.5 * coeffs[1] * vk
+ 0.75 * coeffs[3] * vk ** 3 * i
+ 0.125 * coeffs[5] * vk ** 5 * (10.0 * i ** 2 + 5.0)
)
for i in idx + 1
]
# diagonal wih offset 2, i == j + 2
k = rank - 2
if k > 0:
idx = np.arange(k)
Hamil[idx + 2, idx] = [
factsqrt(i, 2)
* (
0.5 * (uk + coeffs[2] * vk ** 2)
+ 0.25 * coeffs[4] * vk ** 4 * (4.0 * i - 2.0)
+ 15.0 * coeffs[6] * vk ** 6 * (i ** 2 - i + 1.0) / 8.0
)
for i in idx + 2
]
# diagonal wih offset 3, i == j + 3
k = rank - 3
if k > 0:
idx = np.arange(k)
Hamil[idx + 3, idx] = [
factsqrt(i, 3)
* (
coeffs[3] * vk ** 3 / np.sqrt(8.0)
+ np.sqrt(2.0) * coeffs[5] * vk ** 5 * (5.0 * i - 5) / 8.0
)
for i in idx + 3
]
# diagonal wih offset 4, i == j + 4
k = rank - 4
if k > 0:
idx = np.arange(k)
Hamil[idx + 4, idx] = [
factsqrt(i, 4)
* (
0.25e0 * coeffs[4] * vk ** 4
+ 0.125e0 * coeffs[6] * vk ** 6 * (6.0 * i - 9)
)
for i in idx + 4
]
# diagonal wih offset 5, i == j + 5
k = rank - 5
if k > 0:
idx = np.arange(k)
Hamil[idx + 5, idx] = [
np.sqrt(1.0 / 320) * factsqrt(i, 5) * coeffs[5] * vk ** 5 for i in idx + 5
]
# diagonal wih offset 6, i == j + 6
k = rank - 6
if k > 0:
idx = np.arange(k)
Hamil[idx + 6, idx] = [
0.125e0 * factsqrt(i, 6) * coeffs[6] * vk ** 6 for i in idx + 6
]
# make the ful symmetrix matrix
Hamil = Hamil + Hamil.T - np.diag(Hamil.diagonal())
return Hamil
[docs]def anharmonic_frequencies(atoms, T, coeffs, modeinfo):
"""
Calculate the anharmonic frequencies
Parameters
----------
atoms : ase.Atoms
Atoms object
T : float
Temperature in `K`
coeffs : pandas.DataFrame
modeinfo : pandas.DataFrame
"""
MAXITER = 100
QVIB_THRESH = 1.0e-8
FREQ_THRESH = 1.0e-7
au2joule = value("hartree-joule relationship")
invcm2au = 100 * value("inverse meter-hartree relationship")
kT = Boltzmann * T
df = pd.DataFrame(
columns=[
"freq",
"zpve",
"qvib",
"U",
"S",
"converged",
"info",
"rank",
"type",
"d_qvib",
"d_nu",
],
index=modeinfo[modeinfo["vibration"]].index,
dtype=float,
)
for mode, row in coeffs.iterrows():
terminate = False
rank = 4
niter = 0
qvib_last = 0.0
freq_last = 0.0
while not terminate:
# get polynomial coefficients
if row.shape[0] == 7:
pc = row.values[::-1].astype(float)
elif row.shape[0] == 5:
pc = np.append(row.values[::-1], np.zeros(2)).astype(float)
hamil = get_hamiltonian(
rank,
modeinfo.loc[mode, "frequency"] * invcm2au,
modeinfo.loc[mode, "effective_mass"],
pc,
)
w, v = np.linalg.eig(hamil)
w = np.sort(w)
qvib = np.sum(np.exp(-w * au2joule / kT))
anhfreq = (w[1] - w[0]) / invcm2au
zpve = w[0] * au2joule * 1.0e-3 * Avogadro
U, S = get_anh_state_functions(w * au2joule, T)
d_qvib = np.abs(qvib - qvib_last)
d_nu = np.abs(w[0] - freq_last)
terminate = (d_qvib < QVIB_THRESH) & (d_nu < FREQ_THRESH)
if terminate:
if anhfreq < modeinfo.loc[mode, "frequency"]:
anh = (
anhfreq,
zpve,
qvib,
U,
S,
True,
"OK",
rank,
"A",
d_qvib,
d_nu,
)
else:
anh = (
anhfreq,
zpve,
qvib,
U,
S,
True,
"AGTH",
rank,
"A",
d_qvib,
d_nu,
)
else:
rank += 1
qvib_last = qvib
freq_last = w[0]
if niter >= MAXITER:
anh = (
anhfreq,
zpve,
qvib,
U,
S,
False,
"MAXITER",
rank,
"A",
d_qvib,
d_nu,
)
break
niter += 1
df.loc[mode] = anh
df["rank"] = df["rank"].fillna(0).astype(int)
return df
[docs]def merge_vibs(anh6, anh4, harmonic, verbose=False):
"""
Form a DataFrame with the per mode thermochemical
contributions from three separate dataframes with sixth order
polynomial fitted potentia, fourth order fitted potential and
harmonic frequencies.
Parameters
----------
anh6 : pandas.DataFrame
anh4 : pandas.DataFrame
harmonic : pandas.DataFrame
Returns
-------
df : pandas.DataFrame
"""
anh6["order"] = 6
anh4["order"] = 4
if verbose:
print("\n" + " Thermochemistry per mode harmonic ".center(80, "="), end="\n\n")
print_mode_thermo(harmonic)
print("\n" + " Thermochemistry per mode 6th order ".center(80, "="), end="\n\n")
print_mode_thermo(anh6)
print("\n" + " Thermochemistry per mode 4th order ".center(80, "="), end="\n\n")
print_mode_thermo(anh4)
df = pd.DataFrame(columns=anh6.columns, index=anh6.index)
df.update(anh4[anh4["converged"]])
df.update(anh6[anh6["info"] == "OK"])
for col in ["freq", "zpve", "qvib", "U", "S"]:
df[col] = df[col].astype(float)
if df.isnull().any(axis=1).any():
df.fillna(harmonic, inplace=True)
if verbose:
print(
"\n" + " Final data to be used for thermochemistry ".center(80, "="),
end="\n\n",
)
print_mode_thermo(df)
return df
[docs]def harmonic_df(modeinfo, T):
"""
Calculate per mode contributions to the thermodynamic functions in the
harmonic approximation
Parameters
----------
modeinfo : pandas.DataFrame
T : float
Temperature in `K`
Returns
-------
df : pandas.DataFrame
"""
df = pd.DataFrame(
columns=["freq", "zpve", "qvib", "U", "S", "energy", "type"],
index=modeinfo.index,
dtype=float,
)
kT = Boltzmann * T
df["type"] = "H"
df["freq"] = modeinfo["frequency"]
df["energy"] = (
Planck * df["freq"] * 100.0 * value("inverse meter-hertz relationship")
)
df = df[df["freq"] > 0.0]
df["zpve"] = 0.5 * df["energy"] * 1.0e-3 * Avogadro
df["qvib"] = 1.0 / (1.0 - np.exp(-df["energy"] / kT))
df["U"] = (
df["zpve"]
+ 1.0e-3
* gas_constant
* df["energy"]
/ (np.exp(df["energy"] / kT) - 1.0)
/ Boltzmann
)
df["S"] = (
1.0e-3
* gas_constant
* (
df["energy"] / (np.exp(df["energy"] / kT) - 1.0) / kT
- np.log(1.0 - np.exp(-df["energy"] / kT))
)
)
return df
[docs]def get_anh_state_functions(eigenvals, T):
"""
Calculate the internal energy ``U`` and entropy ``S`` for an anharmonic
vibrational mode with eigenvalues ``eigvals`` at temperature ``T`` in
kJ/mol
.. math::
U &= N_{A}\\frac{\sum^{n}_{i=1} \epsilon_{i}\exp(\epsilon_{i}/k_{B}T) }{\sum^{n}_{i=1} \exp(\epsilon_{i}/k_{B}T)}
S &= N_{A}k_{B}\log(\sum^{n}_{i=1} \exp(\epsilon_{i}/k_{B}T))
+ \\frac{N_{A}}{T}\\frac{\sum^{n}_{i=1} \epsilon_{i}\exp(\epsilon_{i}/k_{B}T) }{\sum^{n}_{i=1} \exp(\epsilon_{i}/k_{B}T)}
Parameters
----------
eigenvals : numpy.array
Eigenvalues of the anharmonic 1D Hamiltonian in Joules
T : float
Temperature in `K`
Returns
-------
(U, S) : tuple of floats
Tuple with the internal energy and entropy in kJ/mol
"""
kT = Boltzmann * T
sum1 = np.sum(eigenvals * np.exp(-eigenvals / kT))
sum2 = np.sum(np.exp(-eigenvals / kT))
U = Avogadro * sum1 / sum2
S = Boltzmann * Avogadro * np.log(sum2) + Avogadro * sum1 / (sum2 * T)
# convert J/mol to kJ/mol
return (U * 1.0e-3, S * 1.0e-3)