# Source code for strawberryfields.apps.qchem.utils

# Copyright 2020 Xanadu Quantum Technologies Inc.

# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
r"""
This module contains various utility functions needed to perform quantum chemistry calculations with
Strawberry Fields.
"""
from typing import Tuple

import numpy as np
from scipy.constants import c, h, m_u, pi

[docs]def duschinsky(
Li: np.ndarray, Lf: np.ndarray, ri: np.ndarray, rf: np.ndarray, wf: np.ndarray, m: np.ndarray
) -> Tuple[np.ndarray, np.ndarray]:
r"""Generate the Duschinsky rotation matrix :math:U and displacement vector :math:\delta.

The Duschinsky transformation relates the normal coordinates of the initial and
final states in a vibronic transition, :math:q_i and :math:q_f respectively, as:

.. math:: q_f = U q_i + d,

where :math:U is the Duschinsky rotation matrix and :math:d is a vector giving the
displacement between the equilibrium structures of the two states involved in the vibronic
transition. The normal coordinates of a molecule can be represented in terms of atomic
displacements as:

.. math:: q = L^T \sqrt{m} (r -r_e),

where :math:r_e represents the equilibrium geometry of the molecule, :math:m represents
atomic masses and :math:L is a matrix containing the eigenvectors of the mass-weighted
Hessian. The Duschinsky parameters :math:U and :math:d can be obtained as:

.. math:: U = L_f^T L_i,

.. math:: d = L_f^T \sqrt{m} (r_e^i-r_e^f).

Note that :math:i and :math:f refer to the initial and final states, respectively. The
parameter :math:d is usually represented as a dimensionless parameter :math:\delta as:

.. math:: \delta = l^{-1} d,

where :math:l is a diagonal matrix containing the vibrational frequencies :math:\omega of
the final state:

.. math:: l_{kk} = \left ( \frac{\hbar }{2 \pi \omega_k c} \right )^{1/2},

where :math:\hbar is the reduced Planck constant and :math:c is the speed of light.

The vibrational normal mode matrix for a molecule with :math:M vibrational modes and
:math:N atoms is a :math:3N \times M matrix where :math:M = 3N - 6 for nonlinear molecules and
:math:M = 3N - 5 for linear molecules. The Duschinsky rotation matrix of a molecule is an
:math:M \times M matrix and the Duschinsky displacement vector has :math:M components.

**Example usage:**

>>> Li = np.array([[-0.28933191], [0.0], [0.0], [0.95711104], [0.0], [0.0]])
>>> Lf = np.array([[-0.28933191], [0.0], [0.0], [0.95711104], [0.0], [0.0]])
>>> ri = np.array([-0.0236, 0.0, 0.0, 1.2236, 0.0, 0.0])
>>> rf = np.array([0.0, 0.0, 0.0, 1.4397, 0.0, 0.0])
>>> wf = np.array([1363.2])
>>> m = np.array([11.0093] * 3 + [1.0078] * 3)
>>> U, delta = duschinsky(Li, Lf, ri, rf, wf, m)
>>> U, delta
(array([[0.99977449]]), array([-1.17623073]))

Args:
Li (array): mass-weighted normal modes of the initial electronic state
Lf (array): mass-weighted normal modes of the final electronic state
ri (array): equilibrium molecular geometry of the initial electronic state
rf (array): equilibrium molecular geometry of the final electronic state
wf (array): normal mode frequencies of the final electronic state in units of :math:\mbox{cm}^{-1}
m (array): atomic masses in unified atomic mass units

Returns:
tuple[array, array]: Duschinsky rotation matrix :math:U, Duschinsky displacement vector
:math:\delta
"""
U = Lf.T @ Li

d = Lf.T * m ** 0.5 @ (ri - rf)

l0_inv = np.diag((h / (wf * 100.0 * c)) ** (-0.5) * 2.0 * pi) / 1.0e10 * m_u ** 0.5

delta = np.array(d @ l0_inv)

return U, delta

[docs]def read_gamess(file) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
r"""Reads molecular data from the output file generated by the GAMESS quantum chemistry package
:cite:schmidt1993general.

This function extracts the atomic coordinates (r), atomic masses (m), vibrational frequencies
(w), and normal modes (l) of a molecule from the output file of a vibrational frequency
calculation performed with the GAMESS quantum chemistry package. The output file must contain
the results of a RUNTYP=HESSIAN calculation performed with GAMESS. We recommend checking the
output of this function with the GAMESS results to assure that the GAMESS output file is parsed
correctly.

**Example usage:**

>>> r, m, w, l = read_gamess('../BH_data.out')
>>> r # atomic coordinates
array([[0.0000000, 0.0000000, 0.0000000],
[1.2536039, 0.0000000, 0.0000000]])
>>> m # atomic masses
array([11.00931,  1.00782])
>>> w # vibrational frequencies
array([19.74, 19.73, 0.00, 0.00, 0.00, 2320.32])
>>> l # normal modes
array([[-0.0000000e+00, -7.5322000e-04, -8.7276210e-02,  0.0000000e+00,
8.2280900e-03,  9.5339055e-01],
[-0.0000000e+00, -8.7276210e-02,  7.5322000e-04,  0.0000000e+00,
9.5339055e-01, -8.2280900e-03],
[ 2.8846925e-01, -2.0000000e-08,  2.0000000e-08,  2.8846925e-01,
-2.0000000e-08,  2.0000000e-08],
[ 2.0000000e-08,  2.8846925e-01, -2.0000000e-08,  2.0000000e-08,
2.8846925e-01, -2.0000000e-08],
[-2.0000000e-08,  2.0000000e-08,  2.8846925e-01, -2.0000000e-08,
2.0000000e-08,  2.8846925e-01],
[-8.7279460e-02,  0.0000000e+00,  0.0000000e+00,  9.5342606e-01,
-0.0000000e+00, -0.0000000e+00]])

Args:
file (str): path to the GAMESS output file

Returns:
tuple[array, array, array, array]: atomic coordinates, atomic masses, normal mode
frequencies, normal modes
"""
with open(file, "r") as f:

r = []
m = []
w = []
l = []

for line in f:

if "INPUT CARD> $data" in line or "INPUT CARD>$DATA" in line:
line = [next(f) for _ in range(3)][-1]
while "end" not in line and "END" not in line:
r.append(np.array(line.rstrip().split()[-3:], float))
line = next(f).rstrip()

if "ATOMIC WEIGHTS" in line:
next(f)
for _ in range(len(r)):
m.append(np.array(next(f).rstrip().split()[-1:], float))

if "FREQUENCY" in line:
line = line.rstrip().split()
n_mode = len(line) - 1
w.append(np.array(line[-n_mode:], float))

pass

d = []
for _ in range(len(r) * 3):
l.append(np.array(d, float).T)

if not r:
raise ValueError("No atomic coordinates found in the output file")

if not m:
raise ValueError("No atomic masses found in the output file")

if not w:
raise ValueError("No vibrational frequencies found in the output file")

return (
np.concatenate(r).reshape(len(r), 3),
np.concatenate(m),
np.concatenate(w),
np.concatenate(l),
)