Source code for quantarhei.qm.oscillators.ho
"""Quantarhei package (http://www.github.com/quantarhei)
ho (harmonic oscillator) module
"""
from __future__ import annotations
from typing import Any
import numpy
from ... import COMPLEX, REAL
from ...core.dfunction import DFunction
from ...core.saveable import Saveable
from ...exceptions import BasisError
[docs]
class fcstorage(Saveable):
"""FC factor look-up class
Once Frank-Condon factors for some value of the shift are calculated
they can be stored here, and retrieved when needed again
"""
[docs]
def __init__(self) -> None:
"""Constructor"""
self._shifts: list[float] = []
self._fcs: list[numpy.ndarray] = []
[docs]
def lookup(self, shift: float) -> bool:
"""Returns true if the FC factors for a given shift are available"""
if self._shifts.count(shift) > 0:
return True
return False
[docs]
def index(self, shift: float) -> int:
"""Returns an index of the FC factors with a given shift"""
return self._shifts.index(shift)
[docs]
def add(self, shift: float, fcmatrix: numpy.ndarray) -> None:
"""Adds the matrix of FC factors to the storage"""
self._shifts.append(shift)
self._fcs.append(fcmatrix)
[docs]
def get(self, ii: int) -> numpy.ndarray:
"""Returns a stored FC matrix"""
return self._fcs[ii]
[docs]
class operator_factory(Saveable):
"""Class providing useful operators
Creation and anihilation operators
"""
def __init__(self, N: int = 100) -> None:
# we choose a number of state
# to represent all operators
self.N = N
[docs]
def anihilation_operator(self) -> numpy.ndarray:
N = self.N
aa: numpy.ndarray = numpy.zeros(
(N, N), dtype=REAL
) # matrix N x N full of zeros
for ng in range(N):
for mg in range(N):
if ng == mg - 1:
aa[ng, mg] = numpy.sqrt(numpy.real(mg))
return aa
[docs]
def creation_operator(self) -> numpy.ndarray:
N = self.N
ad: numpy.ndarray = numpy.zeros((N, N), dtype=REAL)
for ng in range(N):
for mg in range(N):
if ng == mg + 1:
ad[ng, mg] = numpy.sqrt(numpy.real(mg + 1))
return ad
[docs]
def shift_operator(self, dd_: Any) -> numpy.ndarray:
"""Calculates the Shift Operator based on the size N_ of the basis
of states and the shift dd_.
The shift operator is defined as
.. math::
D_{\\alpha} = e^{-\\alpha \\frac{\\partial}{\\partial Q}}
where :math:`Q` is the dimensionless coordinate of the Harmonic
oscillator with Hamiltonian
.. math::
H = \\frac{\\hbar\\omega}{2}\\left(P^2 + Q^2\\right).
In this definition, the shift operator acts on a statevector
:math:`\\psi(Q)` in :math:`Q`-representation is such a way that it
shifts it by the value :math:`\\alpha` along the its coordinate
:math:`Q`, i.e.
.. math::
D_{\\alpha}\\psi(Q) = \\psi(Q-\\alpha).
This definition is consistent with the definition of a shift operator
as defined in Wikipedia (look for shift operator). It shifts the
function to the right along the coordite :math:`Q` axis (unlike in
the definition in Wikipedia - this seems to be more natural for
physicists.)
The dimensionless
coordite :math:`Q` and dimensionless momentum :math:`P` are related
to creation and annihilation operators as
.. math::
a = \\frac{1}{\\sqrt{2}}\\left(Q+iP\\right)
.. math::
a^{\\dagger} = \\frac{1}{\\sqrt{2}}\\left(Q-iP\\right)
.. math::
Q = \\frac{1}{\\sqrt{2}}\\left(a + a^{\\dagger}\\right)
.. math::
P = \\frac{1}{i\\sqrt{2}}\\left(a-a^{\\dagger}\\right).
With these definitions we have
.. math::
H = \\hbar\\omega\\left(a^{\\dagger}a + \\frac{1}{2}\\right)
As we have :math:`P=-i\\frac{\\partial}{\\partial Q}`, the shift
operator reads as
.. math::
D_{\\alpha} = e^{-i\\alpha P}
= e^{-\\frac{\\alpha}{\\sqrt{2}}\\left(a-a^{\\dagger}\\right)}
The operator can be generalized for complex values of :math:`\\alpha`
to read
.. math::
D_{\\alpha} =
e^{\\frac{1}{\\sqrt{2}}\\left(\\alpha
a^{\\dagger}-\\alpha^{*}a\\right)},
where :math:`*` represents complex conjugation. This is the definition
implemented in Quantarhei.
The definition differs by the factor of :math:`\\sqrt{2}` from what
is usually used in literature, but for applications in molecular
physics, this definition seems to be more reasonable.
"""
N_ = self.N
aa = self.anihilation_operator()
ad = self.creation_operator()
# construct the Shift Operator
Dd_large: numpy.ndarray = numpy.zeros((N_, N_), dtype=COMPLEX)
Dd_large = (dd_ * ad - numpy.conj(dd_) * aa) / numpy.sqrt(2.0)
# Diagonalize and obtain transformation matrix
A, S = numpy.linalg.eig(Dd_large)
S1 = numpy.linalg.inv(S)
# Exponentiate
Dd_large = numpy.diag(numpy.exp(A))
# Transform back and reduce to the lower number of states
return numpy.dot(S, numpy.dot(Dd_large, S1))
[docs]
def unity_operator(self) -> numpy.ndarray:
ones: numpy.ndarray = numpy.ones(self.N, dtype=REAL)
ret = numpy.diag(ones)
return ret
[docs]
class qrepresentation:
"""Coordinate representation of the HO wavefunctions"""
def __init__(self, qaxis: Any) -> None:
self.qaxis = qaxis
self.ho_eigenfce_generated = False
[docs]
def generate_ho_eigenfunctions(self) -> None:
"""Generated q-representation of HO eigenfunctions"""
self.ho_eigenfce_generated = True
[docs]
def get_ho_eigenfunction(self, N: int) -> None:
"""Returns q-representation of HO eigenfunction"""
if self.ho_eigenfce_generated:
pass
else:
raise BasisError("HO Eigenfunctions must be generated first")
[docs]
def get_ho_ground_state(self) -> DFunction:
"""Returns the ground state wavefunction of the Harmonic oscillator"""
data = numpy.exp(-(self.qaxis.data**2) / 2) / numpy.sqrt(numpy.sqrt(numpy.pi))
psi0 = DFunction(x=self.qaxis, y=data)
return psi0
[docs]
def get_coherent_state(self, alpha: Any) -> DFunction:
"""Returns q-representation of a coherent state with a given alpha"""
ar = numpy.real(alpha)
ai = numpy.imag(alpha)
sq2 = numpy.sqrt(2.0)
q = self.qaxis.data
data = (
numpy.exp(-((q + sq2 * ar) ** 2) / 2)
* numpy.exp(-1j * ar * ai + 1j * sq2 * ai * q)
/ numpy.sqrt(numpy.sqrt(numpy.pi))
)
psi_alpha = DFunction(x=self.qaxis, y=data)
return psi_alpha
[docs]
def get_probability_distribution(self, wfce: Any) -> DFunction:
"""Returns probability distribution for a given wavefunction"""
return DFunction(x=wfce.axis, y=numpy.abs(wfce.data) ** 2)