Source code for pennylane.qchem.vibrational.taylor_ham
# Copyright 2018-2024 Xanadu Quantum Technologies Inc.
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""The functions related to the construction of the Taylor form Hamiltonian."""
import itertools
import numpy as np
from pennylane.bose import BoseSentence, BoseWord, binary_mapping, unary_mapping
# pylint: disable=import-outside-toplevel
def _import_sklearn():
"""Import sklearn."""
try:
import sklearn
except ImportError as Error:
raise ImportError(
"This feature requires sklearn. It can be installed with: pip install scikit-learn."
) from Error
return sklearn
def _remove_harmonic(freqs, onemode_pes):
"""Removes the harmonic part from the PES.
Args:
freqs (list(float)): the harmonic frequencies in atomic units
onemode_pes (TensorLike[float]): one mode PES
Returns:
tuple: A tuple containing the following:
- TensorLike[float] : anharmonic part of the PES
- TensorLike[float] : harmonic part of the PES
"""
nmodes, quad_order = np.shape(onemode_pes)
grid, _ = np.polynomial.hermite.hermgauss(quad_order)
harmonic_pes = np.zeros((nmodes, quad_order))
anh_pes = np.zeros((nmodes, quad_order))
for ii in range(nmodes):
ho_const = freqs[ii] / 2
harmonic_pes[ii, :] = ho_const * (grid**2)
anh_pes[ii, :] = onemode_pes[ii, :] - harmonic_pes[ii, :]
return anh_pes, harmonic_pes
def _fit_onebody(onemode_op, max_deg, min_deg=3):
r"""Fits the one-mode operator to get one-body coefficients.
Args:
onemode_op (TensorLike[float]): one-mode operator
max_deg (int): maximum degree of Taylor form polynomial
min_deg (int): minimum degree of Taylor form polynomial
Returns:
tuple (TensorLike[float], TensorLike[float]):
- the one-body coefficients
- the predicted one-body PES using fitted coefficients
"""
_import_sklearn()
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
if max_deg < min_deg:
raise ValueError(
f"Taylor expansion degree is {max_deg}<{min_deg}, please set max_deg greater than min_deg."
)
nmodes, quad_order = np.shape(onemode_op)
grid, _ = np.polynomial.hermite.hermgauss(quad_order)
coeffs = np.zeros((nmodes, max_deg - min_deg + 1))
predicted_1D = np.zeros_like(onemode_op)
for i1 in range(nmodes):
poly1D = PolynomialFeatures(degree=(min_deg, max_deg), include_bias=False)
poly1D_features = poly1D.fit_transform(grid.reshape(-1, 1))
poly1D_reg_model = LinearRegression()
poly1D_reg_model.fit(poly1D_features, onemode_op[i1, :])
coeffs[i1, :] = poly1D_reg_model.coef_
predicted_1D[i1, :] = poly1D_reg_model.predict(poly1D_features)
return coeffs, predicted_1D
def _twobody_degs(max_deg, min_deg=3):
"""Finds the degree of fit for two-body coefficients.
Args:
max_deg (int): maximum degree of Taylor form polynomial
min_deg (int): minimum degree of Taylor form polynomial
Returns:
list(tuple): A list of tuples `(q1deg, q2deg)` where the sum of the two values is
guaranteed to be between the maximum total degree and minimum degree.
"""
fit_degs = []
for feat_deg in range(min_deg, max_deg + 1):
max_deg = feat_deg - 1
for deg_dist in range(1, max_deg + 1):
q1deg = max_deg - deg_dist + 1
q2deg = deg_dist
fit_degs.append((q1deg, q2deg))
return fit_degs
def _fit_twobody(twomode_op, max_deg, min_deg=3):
r"""Fits the two-mode operator to get two-body coefficients.
Args:
twomode_op (TensorLike[float]): two-mode operator
max_deg (int): maximum degree of Taylor form polynomial
min_deg (int): minimum degree of Taylor form polynomial
Returns:
tuple (TensorLike[float], TensorLike[float]):
- the two-body coefficients
- the predicted two-body PES using fitted coefficients
"""
_import_sklearn()
from sklearn.linear_model import LinearRegression
nmodes, _, quad_order, _ = np.shape(twomode_op)
gauss_grid, _ = np.polynomial.hermite.hermgauss(quad_order)
if max_deg < min_deg:
raise ValueError(
f"Taylor expansion degree is {max_deg}<{min_deg}, please set max_deg greater than min_deg."
)
fit_degs = _twobody_degs(max_deg, min_deg)
num_coeffs = len(fit_degs)
coeffs = np.zeros((nmodes, nmodes, num_coeffs))
predicted_2D = np.zeros_like(twomode_op)
grid_2D = np.array(np.meshgrid(gauss_grid, gauss_grid))
q1, q2 = (grid.flatten() for grid in grid_2D)
idx_2D = np.array(np.meshgrid(range(quad_order), range(quad_order)))
idx1, idx2 = (idx.flatten() for idx in idx_2D)
num_2D = len(q1)
features = np.zeros((num_2D, num_coeffs))
for deg_idx, (q1deg, q2deg) in enumerate(fit_degs):
features[:, deg_idx] = (q1**q1deg) * (q2**q2deg)
for i1 in range(nmodes):
for i2 in range(i1):
Y = twomode_op[i1, i2, idx1, idx2]
poly2D_reg_model = LinearRegression()
poly2D_reg_model.fit(features, Y)
coeffs[i1, i2, :] = poly2D_reg_model.coef_
predicted = poly2D_reg_model.predict(features)
for idx in range(num_2D):
predicted_2D[i1, i2, idx1[idx], idx2[idx]] = predicted[idx]
return coeffs, predicted_2D
def _generate_bin_occupations(max_occ, nbins):
"""
Generate all valid combinations of bin occupations for a given number of bins
and a total maximum occupancy.
Args:
max_occ(int): the maximum total number of items to be distributed across bins
nbins(int): the number of bins to distribute the items into
Returns
list(tuple): A list of tuples, where each tuple represents a valid combination of item counts for the bins.
"""
combinations = list(itertools.product(range(max_occ + 1), repeat=nbins))
valid_combinations = [combo for combo in combinations if sum(combo) == max_occ]
return valid_combinations
def _threebody_degs(max_deg, min_deg=3):
"""Finds the degree of fit for three-body coefficients.
Args:
max_deg (int): maximum degree of Taylor form polynomial
min_deg (int): minimum degree of Taylor form polynomial
Returns:
list(tuple): A list of tuples `(q1deg, q2deg, q3deg)` where the sum of the three values is
guaranteed to be between the maximum total degree and minimum degree.
"""
fit_degs = []
for feat_deg in range(min_deg, max_deg + 1):
max_deg = feat_deg - 3
if max_deg < 0:
continue
possible_occupations = _generate_bin_occupations(max_deg, 3)
for occ in possible_occupations:
q1deg = 1 + occ[0]
q2deg = 1 + occ[1]
q3deg = 1 + occ[2]
fit_degs.append((q1deg, q2deg, q3deg))
return fit_degs
def _fit_threebody(threemode_op, max_deg, min_deg=3):
r"""Fits the three-mode operator to get three-body coefficients.
Args:
threemode_op (TensorLike[float]): threemode operator
max_deg (int): maximum degree of Taylor form polynomial
min_deg (int): minimum degree of Taylor form polynomial
Returns:
tuple (TensorLike[float], TensorLike[float]):
- the three-body coefficients
- the predicted three-body PES using fitted coefficients
"""
_import_sklearn()
from sklearn.linear_model import LinearRegression
nmodes, _, _, quad_order, _, _ = np.shape(threemode_op)
gauss_grid, _ = np.polynomial.hermite.hermgauss(quad_order)
if max_deg < min_deg:
raise ValueError(
f"Taylor expansion degree is {max_deg}<{min_deg}, please set max_deg greater than min_deg."
)
predicted_3D = np.zeros_like(threemode_op)
fit_degs = _threebody_degs(max_deg)
num_coeffs = len(fit_degs)
coeffs = np.zeros((nmodes, nmodes, nmodes, num_coeffs))
grid_3D = np.array(np.meshgrid(gauss_grid, gauss_grid, gauss_grid))
q1, q2, q3 = (grid.flatten() for grid in grid_3D)
idx_3D = np.array(np.meshgrid(range(quad_order), range(quad_order), range(quad_order)))
idx1, idx2, idx3 = (idx.flatten() for idx in idx_3D)
num_3D = len(q1)
features = np.zeros((num_3D, num_coeffs))
for deg_idx, Qs in enumerate(fit_degs):
q1deg, q2deg, q3deg = Qs
features[:, deg_idx] = q1 ** (q1deg) * q2 ** (q2deg) * q3 ** (q3deg)
for i1 in range(nmodes):
for i2 in range(i1):
for i3 in range(i2):
Y = threemode_op[i1, i2, i3, idx1, idx2, idx3]
poly3D_reg_model = LinearRegression()
poly3D_reg_model.fit(features, Y)
coeffs[i1, i2, i3, :] = poly3D_reg_model.coef_
predicted = poly3D_reg_model.predict(features)
for idx in range(num_3D):
predicted_3D[i1, i2, i3, idx1[idx], idx2[idx], idx3[idx]] = predicted[idx]
return coeffs, predicted_3D
[docs]def taylor_coeffs(pes, max_deg=4, min_deg=3):
r"""Compute fitted coefficients for Taylor vibrational Hamiltonian.
The coefficients are defined following Eq. 5 of `arXiv:1703.09313
<https://arxiv.org/abs/1703.09313>`_ as:
.. math::
\Phi_{ijk} = \frac{k_{ijk}}{\sqrt{\omega_i \omega_j \omega_k}}
\quad \text{and} \quad
\Phi_{ijkl} = \frac{k_{ijkl}}{\sqrt{\omega_i \omega_j \omega_k \omega_l}},
where :math:`\Phi_{ijk}` and :math:`\Phi_{ijkl}` are the third- and fourth-order reduced force
constants, respectively, defined in terms of the third- and fourth-order partial derivatives
of the potential energy surface data.
Args:
pes (VibrationalPES): object containing the vibrational potential energy surface data
max_deg (int): maximum degree of Taylor form polynomial
min_deg (int): minimum degree of Taylor form polynomial
Returns:
tuple(TensorLike[float]): the coefficients of the one-body, two-body and three-body terms
**Example**
>>> pes_onemode = np.array([[0.309, 0.115, 0.038, 0.008, 0.000, 0.006, 0.020, 0.041, 0.070]])
>>> pes_twomode = np.zeros((1, 1, 9, 9))
>>> dipole_onemode = np.zeros((1, 9, 3))
>>> gauss_weights = np.array([3.96e-05, 4.94e-03, 8.85e-02,
... 4.33e-01, 7.20e-01, 4.33e-01,
... 8.85e-02, 4.94e-03, 3.96e-05])
>>> grid = np.array([-3.19, -2.27, -1.47, -0.72, 0.0, 0.72, 1.47, 2.27, 3.19])
>>> pes_object = qml.qchem.VibrationalPES(
... freqs=np.array([0.025]),
... grid=grid,
... uloc=np.array([[1.0]]),
... gauss_weights=gauss_weights,
... pes_data=[pes_onemode, pes_twomode],
... dipole_data=[dipole_onemode],
... localized=True,
... dipole_level=1,
... )
>>> one, two = qml.qchem.taylor_coeffs(pes_object, 4, 2)
>>> print(one)
[[-0.00088528 -0.00361425 0.00068143]]
"""
anh_pes, harmonic_pes = _remove_harmonic(pes.freqs, pes.pes_onemode)
coeff_1D, predicted_1D = _fit_onebody(anh_pes, max_deg, min_deg=min_deg)
predicted_1D += harmonic_pes
coeff_arr = [coeff_1D]
predicted_arr = [predicted_1D]
if pes.pes_twomode is not None:
coeff_2D, predicted_2D = _fit_twobody(pes.pes_twomode, max_deg, min_deg=min_deg)
coeff_arr.append(coeff_2D)
predicted_arr.append(predicted_2D)
if pes.pes_threemode is not None:
coeff_3D, predicted_3D = _fit_threebody(pes.pes_threemode, max_deg, min_deg=min_deg)
coeff_arr.append(coeff_3D)
predicted_arr.append(predicted_3D)
return coeff_arr
[docs]def taylor_dipole_coeffs(pes, max_deg=4, min_deg=1):
r"""Compute fitted coefficients for the Taylor dipole operator.
Args:
pes (VibrationalPES): object containing the vibrational potential energy surface data
max_deg (int): maximum degree of Taylor form polynomial
min_deg (int): minimum degree of Taylor form polynomial
Returns:
tuple: a tuple containing:
- list(array(floats)): coefficients for x-displacements
- list(array(floats)): coefficients for y-displacements
- list(array(floats)): coefficients for z-displacements
**Example**
>>> freqs = np.array([0.01885397])
>>> grid, weights = np.polynomial.hermite.hermgauss(9)
>>> pes_onebody = np.array([[0.05235573, 0.03093067, 0.01501878, 0.00420778, 0.0,
... 0.00584504, 0.02881817, 0.08483433, 0.22025702]])
>>> pes_twobody = None
>>> dipole_onebody = np.array([[[-1.92201700e-16, 1.45397041e-16, -1.40451549e-01],
... [-1.51005108e-16, 9.53185441e-17, -1.03377032e-01],
... [-1.22793018e-16, 7.22781963e-17, -6.92825934e-02],
... [-1.96537436e-16, -5.86686504e-19, -3.52245369e-02],
... [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
... [ 5.24758835e-17, -1.40650833e-16, 3.69955543e-02],
... [-4.52407941e-17, 1.38406311e-16, 7.60888733e-02],
... [-4.63820104e-16, 5.42928787e-17, 1.17726042e-01],
... [ 1.19224372e-16, 9.12491386e-17, 1.64013197e-01]]])
>>> vib_obj = qml.qchem.VibrationalPES(
... freqs=freqs,
... grid=grid,
... gauss_weights=weights,
... uloc=None,
... pes_data=[pes_onebody, pes_twobody],
... dipole_data=[dipole_onebody],
... localized=False
... )
>>> x, y, z = qml.qchem.taylor_dipole_coeffs(vib_obj, 4, 2)
>>> print(z)
[array([[ 1.64124324e-03, 5.39120159e-03, -4.80053702e-05]])]
"""
coeffs_x_1D, predicted_x_1D = _fit_onebody(
pes.dipole_onemode[:, :, 0], max_deg, min_deg=min_deg
)
coeffs_x_arr = [coeffs_x_1D]
predicted_x_arr = [predicted_x_1D]
coeffs_y_1D, predicted_y_1D = _fit_onebody(
pes.dipole_onemode[:, :, 1], max_deg, min_deg=min_deg
)
coeffs_y_arr = [coeffs_y_1D]
predicted_y_arr = [predicted_y_1D]
coeffs_z_1D, predicted_z_1D = _fit_onebody(
pes.dipole_onemode[:, :, 2], max_deg, min_deg=min_deg
)
coeffs_z_arr = [coeffs_z_1D]
predicted_z_arr = [predicted_z_1D]
if pes.dipole_twomode is not None:
coeffs_x_2D, predicted_x_2D = _fit_twobody(
pes.dipole_twomode[:, :, :, :, 0], max_deg, min_deg=min_deg
)
coeffs_x_arr.append(coeffs_x_2D)
predicted_x_arr.append(predicted_x_2D)
coeffs_y_2D, predicted_y_2D = _fit_twobody(
pes.dipole_twomode[:, :, :, :, 1], max_deg, min_deg=min_deg
)
coeffs_y_arr.append(coeffs_y_2D)
predicted_y_arr.append(predicted_y_2D)
coeffs_z_2D, predicted_z_2D = _fit_twobody(
pes.dipole_twomode[:, :, :, :, 2], max_deg, min_deg=min_deg
)
coeffs_z_arr.append(coeffs_z_2D)
predicted_z_arr.append(predicted_z_2D)
if pes.dipole_threemode is not None:
coeffs_x_3D, predicted_x_3D = _fit_threebody(
pes.dipole_threemode[:, :, :, :, :, :, 0], max_deg, min_deg=min_deg
)
coeffs_x_arr.append(coeffs_x_3D)
predicted_x_arr.append(predicted_x_3D)
coeffs_y_3D, predicted_y_3D = _fit_threebody(
pes.dipole_threemode[:, :, :, :, :, :, 1], max_deg, min_deg=min_deg
)
coeffs_y_arr.append(coeffs_y_3D)
predicted_y_arr.append(predicted_y_3D)
coeffs_z_3D, predicted_z_3D = _fit_threebody(
pes.dipole_threemode[:, :, :, :, :, :, 2], max_deg, min_deg=min_deg
)
coeffs_z_arr.append(coeffs_z_3D)
predicted_z_arr.append(predicted_z_3D)
return coeffs_x_arr, coeffs_y_arr, coeffs_z_arr
def _position_to_boson(index, op):
"""Convert position operator `p` or `q` into respective bosonic operator. The conversion is
described in `Eq. 6 and 7 <https://arxiv.org/pdf/1703.09313>`_.
Args:
index (int): the index of the operator
op (str): the position operator, either ``"p"`` or ``"q"``
Returns:
BoseSentence: bosonic form of the given position operator
"""
factor = 1j / np.sqrt(2) if op == "p" else 1 / np.sqrt(2)
bop = factor * BoseWord({(0, index): "-"})
bdag = factor * BoseWord({(0, index): "+"})
return bdag - bop if op == "p" else bdag + bop
def _taylor_anharmonic(taylor_coeffs_array, start_deg=2):
"""Build anharmonic term of Taylor form bosonic observable from provided coefficients described
in `Eq. 10 <https://arxiv.org/pdf/1703.09313>`_.
Args:
taylor_coeffs_array (list(float)): the coeffs of the Taylor expansion
start_deg (int): the starting degree
Returns:
pennylane.bose.BoseSentence: anharmonic part of the Taylor hamiltonian for given coeffs
"""
num_coups = len(taylor_coeffs_array)
taylor_1D = taylor_coeffs_array[0]
num_modes, num_1D_coeffs = np.shape(taylor_1D)
taylor_deg = num_1D_coeffs + start_deg - 1
ordered_dict = BoseSentence({})
# One-mode expansion
for mode in range(num_modes):
bosonized_qm = _position_to_boson(mode, "q")
for deg_i in range(start_deg, taylor_deg + 1):
coeff = taylor_1D[mode, deg_i - start_deg]
qpow = bosonized_qm**deg_i
ordered_dict += (coeff * qpow).normal_order()
# Two-mode expansion
if num_coups > 1:
taylor_2D = taylor_coeffs_array[1]
degs_2d = _twobody_degs(taylor_deg, min_deg=start_deg)
for m1 in range(num_modes):
bosonized_qm1 = _position_to_boson(m1, "q")
for m2 in range(m1):
bosonized_qm2 = _position_to_boson(m2, "q")
for deg_idx, Qs in enumerate(degs_2d):
q1deg, q2deg = Qs[:2]
coeff = taylor_2D[m1, m2, deg_idx]
bosonized_qm1_pow = bosonized_qm1**q1deg
bosonized_qm2_pow = bosonized_qm2**q2deg
ordered_dict += (coeff * bosonized_qm1_pow * bosonized_qm2_pow).normal_order()
# Three-mode expansion
if num_coups > 2:
degs_3d = _threebody_degs(taylor_deg, min_deg=start_deg)
taylor_3D = taylor_coeffs_array[2]
for m1 in range(num_modes):
bosonized_qm1 = _position_to_boson(m1, "q")
for m2 in range(m1):
bosonized_qm2 = _position_to_boson(m2, "q")
for m3 in range(m2):
bosonized_qm3 = _position_to_boson(m3, "q")
for deg_idx, Qs in enumerate(degs_3d):
q1deg, q2deg, q3deg = Qs[:3]
coeff = taylor_3D[m1, m2, m3, deg_idx]
bosonized_qm1_pow = bosonized_qm1**q1deg
bosonized_qm2_pow = bosonized_qm2**q2deg
bosonized_qm3_pow = bosonized_qm3**q3deg
ordered_dict += (
coeff * bosonized_qm1_pow * bosonized_qm2_pow * bosonized_qm3_pow
).normal_order()
return BoseSentence(ordered_dict).normal_order()
def _taylor_kinetic(taylor_coeffs_array, freqs, is_local=True, uloc=None):
"""Build kinetic term of Taylor form bosonic observable from provided coefficients
Args:
taylor_coeffs_array (list(float)): the coeffs of the Taylor expansion
freqs (list(float)): the harmonic frequencies in atomic units
is_local (bool): Flag whether the vibrational modes are localized. Default is ``True``.
uloc (list(list(float))): localization matrix indicating the relationship between original
and localized modes
Returns:
pennylane.bose.BoseSentence: kinetic term of the Taylor hamiltonian for given coeffs
"""
taylor_1D = taylor_coeffs_array[0]
num_modes, _ = np.shape(taylor_1D)
if is_local:
alphas_arr = np.einsum("ij,ik,j,k->jk", uloc, uloc, np.sqrt(freqs), np.sqrt(freqs))
else:
alphas_arr = np.zeros((num_modes, num_modes))
for m in range(num_modes):
alphas_arr[m, m] = freqs[m]
kin_ham = BoseSentence({})
for m1 in range(num_modes):
pm1 = _position_to_boson(m1, "p")
for m2 in range(num_modes):
pm2 = _position_to_boson(m2, "p")
kin_ham += (0.5 * alphas_arr[m1, m2]) * (pm1 * pm2).normal_order()
return kin_ham.normal_order()
def _taylor_harmonic(taylor_coeffs_array, freqs):
"""Build harmonic term of Taylor form bosonic observable from provided coefficients, see first
term of `Eq. 4 and Eq. 7 <https://arxiv.org/pdf/1703.09313>`_.
Args:
taylor_coeffs_array (list(float)): the coeffs of the Taylor expansion
freqs (list(float)): the harmonic frequencies in atomic units
Returns:
pennylane.bose.BoseSentence: harmonic term of the Taylor hamiltonian for given coeffs
"""
taylor_1D = taylor_coeffs_array[0]
num_modes, _ = np.shape(taylor_1D)
harm_pot = BoseSentence({})
for mode in range(num_modes):
bosonized_qm2 = (
_position_to_boson(mode, "q") * _position_to_boson(mode, "q")
).normal_order()
harm_pot += bosonized_qm2 * freqs[mode] * 0.5
return harm_pot.normal_order()
[docs]def taylor_bosonic(coeffs, freqs, is_local=True, uloc=None):
"""Return Taylor bosonic vibrational Hamiltonian.
The construction of the Hamiltonian is based on Eqs. 4-7 of `arXiv:1703.09313 <https://arxiv.org/abs/1703.09313>`_.
Args:
coeffs (list(float)): the coefficients of the Hamiltonian
freqs (list(float)): the harmonic frequencies in atomic units
is_local (bool): Flag whether the vibrational modes are localized. Default is ``True``.
uloc (list(list(float))): localization matrix indicating the relationship between original
and localized modes
Returns:
pennylane.bose.BoseSentence: Taylor bosonic hamiltonian
**Example**
>>> one_mode = np.array([[-0.00088528, -0.00361425, 0.00068143]])
>>> two_mode = np.array([[[0., 0., 0., 0., 0., 0.]]])
>>> freqs = np.array([0.025])
>>> uloc = np.array([[1.0]])
>>> ham = qml.qchem.taylor_bosonic(coeffs=[one_mode, two_mode], freqs=freqs, uloc=uloc)
>>> print(ham)
-0.0012778303419517393 * b⁺(0) b⁺(0) b⁺(0)
+ -0.0038334910258552178 * b⁺(0) b⁺(0) b(0)
+ -0.0038334910258552178 * b⁺(0)
+ -0.0038334910258552178 * b⁺(0) b(0) b(0)
+ -0.0038334910258552178 * b(0)
+ -0.0012778303419517393 * b(0) b(0) b(0)
+ (0.0005795050000000001+0j) * b⁺(0) b⁺(0)
+ (0.026159009999999996+0j) * b⁺(0) b(0)
+ (0.012568432499999997+0j) * I
+ (0.0005795050000000001+0j) * b(0) b(0)
+ 0.00017035749999999995 * b⁺(0) b⁺(0) b⁺(0) b⁺(0)
+ 0.0006814299999999998 * b⁺(0) b⁺(0) b⁺(0) b(0)
+ 0.0010221449999999997 * b⁺(0) b⁺(0) b(0) b(0)
+ 0.0006814299999999998 * b⁺(0) b(0) b(0) b(0)
+ 0.00017035749999999995 * b(0) b(0) b(0) b(0)
"""
if is_local:
start_deg = 2
else:
start_deg = 3
harm_pot = _taylor_harmonic(coeffs, freqs)
ham = _taylor_anharmonic(coeffs, start_deg) + harm_pot
kin_ham = _taylor_kinetic(coeffs, freqs, is_local, uloc)
ham += kin_ham
return ham.normal_order()
# pylint: disable=too-many-positional-arguments, too-many-arguments
[docs]def taylor_hamiltonian(
pes, max_deg=4, min_deg=3, mapping="binary", n_states=2, wire_map=None, tol=1e-12
):
"""Return Taylor vibrational Hamiltonian.
The construction of the Hamiltonian is based on Eqs. 4-7 of `arXiv:1703.09313 <https://arxiv.org/abs/1703.09313>`_.
The Hamiltonian is then converted to a qubit operator with a selected ``mapping`` method.
Args:
pes (VibrationalPES): object containing the vibrational potential energy surface data
max_deg (int): maximum degree of Taylor form polynomial
min_deg (int): minimum degree of Taylor form polynomial
mapping (str): Method used to map to qubit basis. Input values can be ``"binary"``
or ``"unary"``. Default is ``"binary"``.
n_states(int): maximum number of allowed bosonic states
wire_map (dict): A dictionary defining how to map the states of the Bose operator to qubit
wires. If ``None``, integers used to label the bosonic states will be used as wire labels.
Defaults to ``None``.
tol (float): tolerance for discarding the imaginary part of the coefficients
Returns:
Operator: the Taylor Hamiltonian
**Example**
>>> pes_onemode = np.array([[0.309, 0.115, 0.038, 0.008, 0.000, 0.006, 0.020, 0.041, 0.070]])
>>> pes_twomode = np.zeros((1, 1, 9, 9))
>>> dipole_onemode = np.zeros((1, 9, 3))
>>> gauss_weights = np.array([3.96e-05, 4.94e-03, 8.85e-02,
... 4.33e-01, 7.20e-01, 4.33e-01,
... 8.85e-02, 4.94e-03, 3.96e-05])
>>> grid = np.array([-3.19, -2.27, -1.47, -0.72, 0.0, 0.72, 1.47, 2.27, 3.19])
>>> pes_object = qml.qchem.VibrationalPES(
... freqs=np.array([0.025]),
... grid=grid,
... uloc=np.array([[1.0]]),
... gauss_weights=gauss_weights,
... pes_data=[pes_onemode, pes_twomode],
... dipole_data=[dipole_onemode],
... localized=True,
... dipole_level=1,
... )
>>> qml.qchem.taylor_hamiltonian(pes_object, 4, 2)
(
-0.003833496032473659 * X(0)
+ (0.0256479442871582+0j) * I(0)
+ (-0.013079509779221888+0j) * Z(0)
)
"""
mapping.lower().strip()
if mapping not in ["binary", "unary"]:
raise ValueError(
f"Specified mapping {mapping}, is not found. Please use either 'binary' or 'unary' mapping."
)
coeffs_arr = taylor_coeffs(pes, max_deg, min_deg)
bose_op = taylor_bosonic(coeffs_arr, pes.freqs, is_local=pes.localized, uloc=pes.uloc)
if mapping == "binary":
ham = binary_mapping(bose_operator=bose_op, n_states=n_states, wire_map=wire_map, tol=tol)
elif mapping == "unary":
ham = unary_mapping(bose_operator=bose_op, n_states=n_states, wire_map=wire_map, tol=tol)
return ham
_modules/pennylane/qchem/vibrational/taylor_ham
Download Python script
Download Notebook
View on GitHub