Source code for pennylane.qchem.vibrational.vscf

# 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.
"""This module contains functions to perform VSCF calculation."""

import numpy as np

# pylint: disable=too-many-arguments, too-many-positional-arguments


def _build_fock(mode, active_ham_terms, active_terms, modals, h_mat, mode_rot):
    r"""Builds the fock matrix for each vibrational mode.

    Args:
        mode (int): index of the active mode for which the Fock matrix is being calculated
        active_ham_terms (dict): dictionary containing the active Hamiltonian terms
        active_terms (list): list containing the active Hamiltonian terms for the given mode
        modals (list[int]): list containing the maximum number of modals to consider for each mode
        h_mat (array[array[float]]): the Hamiltonian matrix
        mode_rot (TensorLike[float]): rotation matrix for the given vibrational mode

    Returns:
        array[array[float]]: fock matrix for the given mode
    """

    fock_matrix = np.zeros((modals[mode], modals[mode]))
    for term in active_terms:
        ham_term = active_ham_terms[term]
        h_val, modal_indices, excitation_indices = ham_term[0], list(ham_term[1]), list(ham_term[2])

        modal_idx = modal_indices.index(mode)
        i = excitation_indices[modal_idx]
        j = excitation_indices[modal_idx + len(modal_indices)]
        if i < modals[mode] and j < modals[mode]:
            h_val = h_val * np.prod(
                [h_mat[term, modal] for modal in modal_indices if modal != mode]
            )
            fock_matrix[i, j] += h_val

    fock_matrix = mode_rot.T @ fock_matrix @ mode_rot
    return fock_matrix


def _update_h(h_mat, mode, active_ham_terms, mode_rot, active_terms, modals):
    r"""Updates value of Hamiltonian matrix for a mode after other modes have been rotated.

    Args:
        h_mat (array[array[float]]): the Hamiltonian matrix
        mode (int): index of the active mode for which the Hamiltonian matrix is being updated
        active_ham_terms (dict): dictionary containing the active Hamiltonian terms
        mode_rot (array[array[float]]): rotation matrix for a given mode
        active_terms (list): list containing active Hamiltonian terms for the given mode
        modals (list[int]): list containing the maximum number of modals to consider for each mode

    Returns:
        array[array[float]]: the updated Hamiltonian matrix

    """
    u_coeffs = mode_rot[:, 0]

    for t in active_terms:
        ham_term = active_ham_terms[t]
        modal_indices, excitation_indices = ham_term[1], ham_term[2]

        modal_idx = modal_indices.index(mode)
        i = excitation_indices[modal_idx]
        j = excitation_indices[modal_idx + len(modal_indices)]

        if i < modals[mode] and j < modals[mode]:
            h_mat[t, mode] = u_coeffs[i] * u_coeffs[j]

    return h_mat


def _find_active_terms(h_integrals, modals, cutoff):
    r"""Identifies the active terms in the Hamiltonian, following the equations 20-22
    in `J. Chem. Theory Comput. 2010, 6, 235–248 <https://pubs.acs.org/doi/10.1021/ct9004454>`_.
    The equation suggests that if mode m is not contained in a Hamiltonian term, it evaluates to zero.

    Args:
        h_integrals (list[TensorLike[float]]): list of n-mode expansion of Hamiltonian integrals
        modals (list[int]): list containing the maximum number of modals to consider for each mode
        cutoff (float): threshold value for including matrix elements into operator

    Returns:
        tuple: A tuple containing the following:
            - dict: dictionary containing the active Hamiltonian terms
            - dict: dictionary containing the active Hamiltonian terms for all vibrational modes
            - int: total number of active terms in the Hamiltonian

    """
    active_ham_terms = {}
    active_num = 0
    nmodes = np.shape(h_integrals[0])[0]

    for n, ham_n in enumerate(h_integrals):
        for idx, h_val in np.ndenumerate(ham_n):
            if np.abs(h_val) < cutoff:
                continue

            modal_indices = idx[: n + 1]

            if all(modal_indices[i] > modal_indices[i + 1] for i in range(n)):
                excitation_indices = idx[n + 1 :]
                exc_in_modals = True
                for m_idx, m in enumerate(modal_indices):
                    i = excitation_indices[m_idx]
                    j = excitation_indices[m_idx + len(modal_indices)]
                    if i >= modals[m] or j >= modals[m]:
                        exc_in_modals = False
                        break
                if exc_in_modals:
                    active_ham_terms[active_num] = [h_val, modal_indices, excitation_indices]
                    active_num += 1

    active_mode_terms = {
        m: [t for t in range(active_num) if m in active_ham_terms[t][1]] for m in range(nmodes)
    }

    return active_ham_terms, active_mode_terms, active_num


def _fock_energy(h_mat, active_ham_terms, active_mode_terms, modals, mode_rots):
    r"""Calculates vibrational energy.

    Args:
        h_mat (array[array[float]]): the Hamiltonian matrix
        active_ham_terms (dict): dictionary containing the active Hamiltonian terms
        active_mode_terms (dict): list containing the active Hamiltonian terms for all vibrational modes
        modals (list[int]): list containing the maximum number of modals to consider for each mode
        mode_rots (List[TensorLike[float]]): list of rotation matrices for each vibrational mode

    Returns:
        float: vibrational energy

    """
    nmodes = h_mat.shape[1]
    e0s = np.zeros(nmodes)
    for mode in range(nmodes):
        fock_mat = _build_fock(
            mode, active_ham_terms, active_mode_terms[mode], modals, h_mat, mode_rots[mode]
        )
        e0s[mode] = fock_mat[0, 0]

    return np.sum(e0s)


def _vscf(h_integrals, modals, cutoff, tol=1e-8, max_iters=10000):
    r"""Performs the VSCF calculation.

    Args:
        h_integrals (list(TensorLike[float])): list containing Hamiltonian integral matrices
        modals (list(int)): list containing the maximum number of modals to consider for each mode
        cutoff (float): threshold value for including matrix elements into operator.
        tol (float): convergence tolerance for vscf calculation
        max_iters (int): maximum number of iterations for vscf to converge

    Returns:
        tuple: A tuple of the following:
            - float: vscf energy
            - list[TensorLike[float]]: list of rotation matrices for all vibrational modes

    """

    nmodes = np.shape(h_integrals[0])[0]

    active_ham_terms, active_mode_terms, active_num = _find_active_terms(
        h_integrals, modals, cutoff
    )
    mode_rots = [np.eye(modals[mode]) for mode in range(nmodes)]
    h_mat = np.zeros((active_num, nmodes))
    for mode in range(nmodes):
        h_mat = _update_h(
            h_mat, mode, active_ham_terms, mode_rots[mode], active_mode_terms[mode], modals
        )

    e0 = _fock_energy(h_mat, active_ham_terms, active_mode_terms, modals, mode_rots)

    enew = e0 + 2 * tol
    for curr_iter in range(max_iters):
        if curr_iter != 0:
            e0 = enew

        for mode in range(nmodes):
            fock = _build_fock(
                mode, active_ham_terms, active_mode_terms[mode], modals, h_mat, mode_rots[mode]
            )
            _, eigvec = np.linalg.eigh(fock)
            mode_rots[mode] = mode_rots[mode] @ eigvec
            h_mat = _update_h(
                h_mat, mode, active_ham_terms, mode_rots[mode], active_mode_terms[mode], modals
            )
            fock = np.transpose(eigvec) @ fock @ eigvec

        enew = _fock_energy(h_mat, active_ham_terms, active_mode_terms, modals, mode_rots)

        if np.abs(enew - e0) <= tol:
            break

    return enew, mode_rots


def _rotate_one_body(h1, nmodes, mode_rots, modals):
    r"""Rotates one body integrals.

    Args:
        h1 (TensorLike[float]): one-body integrals
        nmodes (int): number of vibrational modes
        mode_rots (list[TensorLike[float]]): list of rotation matrices for all vibrational modes
        modals (list[int]): list containing the maximum number of modals to consider for each mode

    Returns:
        TensorLike[float]: rotated one-body integrals

    """
    imax = np.max(modals)
    h1_rot = np.zeros((nmodes, imax, imax))
    for m in range(nmodes):
        h1_rot[m, : modals[m], : modals[m]] = np.einsum(
            "ij,ia,jb->ab", h1[m, :, :], mode_rots[m][:, : modals[m]], mode_rots[m][:, : modals[m]]
        )

    return h1_rot


def _rotate_two_body(h2, nmodes, mode_rots, modals):
    r"""Rotates two body integrals.

    Args:
        h2 (TensorLike[float]): two-body integrals
        nmodes (int): number of vibrational modes
        mode_rots (list[TensorLike[float]]): list of rotation matrices for all vibrational modes
        modals (list[int]): list containing the maximum number of modals to consider for each mode

    Returns:
        TensorLike[float]: rotated two-body integrals

    """
    imax = np.max(modals)

    U_mats = [mode_rots[m] for m in range(nmodes)]
    h2_rot = np.zeros((nmodes, nmodes, imax, imax, imax, imax))
    for m1 in range(nmodes):
        for m2 in range(nmodes):
            h2_rot[m1, m2, : modals[m1], : modals[m2], : modals[m1], : modals[m2]] = np.einsum(
                "ijkl,ia,jb,kc,ld->abcd",
                h2[m1, m2, :, :, :, :],
                U_mats[m1][:, : modals[m1]],
                U_mats[m2][:, : modals[m2]],
                U_mats[m1][:, : modals[m1]],
                U_mats[m2][:, : modals[m2]],
            )

    return h2_rot


def _rotate_three_body(h3, nmodes, mode_rots, modals):
    r"""Rotates three body integrals.

    Args:
        h3 (TensorLike[float]): three-body integrals
        nmodes (int): number of vibrational modes
        mode_rots (list[TensorLike[float]]): list of rotation matrices for all vibrational modes
        modals (list[int]): list containing the maximum number of modals to consider for each mode

    Returns:
        TensorLike[float]: rotated three-body integrals

    """
    imax = np.max(modals)

    h3_rot = np.zeros((nmodes, nmodes, nmodes, imax, imax, imax, imax, imax, imax))
    for m1 in range(nmodes):
        for m2 in range(nmodes):
            for m3 in range(nmodes):
                h3_rot[
                    m1,
                    m2,
                    m3,
                    : modals[m1],
                    : modals[m2],
                    : modals[m3],
                    : modals[m1],
                    : modals[m2],
                    : modals[m3],
                ] = np.einsum(
                    "ijklmn,ia,jb,kc,ld,me,nf->abcdef",
                    h3[m1, m2, m3, :, :, :, :, :, :],
                    mode_rots[m1][:, : modals[m1]],
                    mode_rots[m2][:, : modals[m2]],
                    mode_rots[m3][:, : modals[m3]],
                    mode_rots[m1][:, : modals[m1]],
                    mode_rots[m2][:, : modals[m2]],
                    mode_rots[m3][:, : modals[m3]],
                )

    return h3_rot


def _rotate_dipole(d_integrals, mode_rots, modals):
    r"""Generates VSCF rotated dipole.

    Args:
        d_integrals (list[TensorLike[float]]): list of n-mode expansion of dipole integrals
        mode_rots (list[TensorLike[float]]): list of rotation matrices for all vibrational modes
        modals (list[int]): list containing the maximum number of modals to consider for each mode

    Returns:
        tuple(TensorLike[float]): a tuple of rotated dipole integrals

    """
    n = len(d_integrals)

    nmodes = np.shape(d_integrals[0])[0]
    imax = np.max(modals)
    d1_rot = np.zeros((3, nmodes, imax, imax))

    for alpha in range(3):
        d1_rot[alpha, ::] = _rotate_one_body(
            d_integrals[0][alpha, ::], nmodes, mode_rots, modals=modals
        )
    dip_data = [d1_rot]

    if n > 1:
        d2_rot = np.zeros((3, nmodes, nmodes, imax, imax, imax, imax))
        for alpha in range(3):
            d2_rot[alpha, ::] = _rotate_two_body(
                d_integrals[1][alpha, ::], nmodes, mode_rots, modals=modals
            )
        dip_data = [d1_rot, d2_rot]

    if n > 2:
        d3_rot = np.zeros((3, nmodes, nmodes, nmodes, imax, imax, imax, imax, imax, imax))
        for alpha in range(3):
            d3_rot[alpha, ::] = _rotate_three_body(
                d_integrals[2][alpha, ::], nmodes, mode_rots, modals=modals
            )
        dip_data = [d1_rot, d2_rot, d3_rot]

    return dip_data


def _rotate_hamiltonian(h_integrals, mode_rots, modals):
    r"""Generates VSCF rotated Hamiltonian.

    Args:
        h_integrals (list[TensorLike[float]]): list of n-mode expansion of Hamiltonian integrals
        mode_rots (list[TensorLike[float]]): list of rotation matrices for all vibrational modes
        modals (list[int]): list containing the maximum number of modals to consider for each mode

    Returns:
        tuple(TensorLike[float]): tuple of rotated Hamiltonian integrals

    """

    n = len(h_integrals)
    nmodes = np.shape(h_integrals[0])[0]

    h1_rot = _rotate_one_body(h_integrals[0], nmodes, mode_rots, modals)
    h_data = [h1_rot]

    if n > 1:
        h2_rot = _rotate_two_body(h_integrals[1], nmodes, mode_rots, modals)
        h_data = [h1_rot, h2_rot]

    if n > 2:
        h3_rot = _rotate_three_body(h_integrals[2], nmodes, mode_rots, modals)
        h_data = [h1_rot, h2_rot, h3_rot]

    return h_data


def vscf_integrals(h_integrals, d_integrals=None, modals=None, cutoff=None, cutoff_ratio=1e-6):
    r"""Generates vibrational self-consistent field rotated integrals.

    This function converts the Christiansen vibrational Hamiltonian integrals obtained in the harmonic
    oscillator basis to integrals in the vibrational self-consistent field (VSCF) basis.
    The implementation is based on the method described in
    `J. Chem. Theory Comput. 2010, 6, 235–248 <https://pubs.acs.org/doi/10.1021/ct9004454>`_.

    Args:
        h_integrals (list[TensorLike[float]]): List of Hamiltonian integrals for up to 3 coupled vibrational modes.
            Look at the Usage Details for more information.
        d_integrals (list[TensorLike[float]]): List of dipole integrals for up to 3 coupled vibrational modes.
            Look at the Usage Details for more information.
        modals (list[int]): list containing the maximum number of modals to consider for each vibrational mode.
            Default value is the maximum number of modals.
        cutoff (float): threshold value for including matrix elements into operator
        cutoff_ratio (float): ratio for discarding elements with respect to biggest element in the integrals.
            Default value is ``1e-6``.

    Returns:
        tuple: a tuple containing:
            - list[TensorLike[float]]: List of Hamiltonian integrals in VSCF basis for up to 3 coupled vibrational modes.
            - list[TensorLike[float]]: List of dipole integrals in VSCF basis for up to 3 coupled vibrational modes.

        ``None`` is returned if ``d_integrals`` is ``None``.

    **Example**

    >>> h1 = np.array([[[0.00968289, 0.00233724, 0.0007408,  0.00199125],
    ...                 [0.00233724, 0.02958449, 0.00675431, 0.0021936],
    ...                 [0.0007408,  0.00675431, 0.0506012,  0.01280986],
    ...                 [0.00199125, 0.0021936,  0.01280986, 0.07282307]]])
    >>> qml.qchem.vscf_integrals(h_integrals=[h1], modals=[4,4,4])
    ([array([[[ 9.36124041e-03, -4.20128342e-19,  3.25260652e-19,
            1.08420217e-18],
            [-9.21571847e-19,  2.77803512e-02, -3.46944695e-18,
            5.63785130e-18],
            [-3.25260652e-19, -8.67361738e-19,  4.63297357e-02,
            -1.04083409e-17],
            [ 1.30104261e-18,  5.20417043e-18, -1.38777878e-17,
            7.92203227e-02]]])],
    None)


    .. details::
        :title: Usage Details

        The ``h_integral`` tensor must have one of these dimensions:

        - 1-mode coupled integrals: `(n, m)`
        - 2-mode coupled integrals: `(n, n, m, m, m, m)`
        - 3-mode coupled integrals: `(n, n, n, m, m, m, m, m, m)`

        where ``n`` is the number of vibrational modes in the molecule and ``m`` represents the number
        of modals.

        The ``d_integral`` tensor must have one of these dimensions:

        - 1-mode coupled integrals: `(3, n, m)`
        - 2-mode coupled integrals: `(3, n, n, m, m, m, m)`
        - 3-mode coupled integrals: `(3, n, n, n, m, m, m, m, m, m)`

        where ``n`` is the number of vibrational modes in the molecule, ``m`` represents the number
        of modals and the first axis represents the ``x, y, z`` component of the dipole. Default is ``None``.

    """

    if len(h_integrals) > 3:
        raise ValueError(
            f"Building n-mode Hamiltonian is not implemented for n equal to {len(h_integrals)}."
        )

    if d_integrals is not None:
        if len(d_integrals) > 3:
            raise ValueError(
                f"Building n-mode dipole is not implemented for n equal to {len(d_integrals)}."
            )

    nmodes = np.shape(h_integrals[0])[0]

    imax = np.shape(h_integrals[0])[1]
    max_modals = nmodes * [imax]
    if modals is None:
        modals = max_modals
    else:
        if np.max(modals) > imax:
            raise ValueError(
                "Number of maximum modals cannot be greater than the modals for unrotated integrals."
            )
        imax = np.max(modals)

    if cutoff is None:
        max_val = np.max([np.max(np.abs(H)) for H in h_integrals])
        cutoff = max_val * cutoff_ratio

    _, mode_rots = _vscf(h_integrals, modals=max_modals, cutoff=cutoff)

    h_data = _rotate_hamiltonian(h_integrals, mode_rots, modals)

    if d_integrals is not None:
        dip_data = _rotate_dipole(d_integrals, mode_rots, modals)
        return h_data, dip_data

    return h_data, None

_modules/pennylane/qchem/vibrational/vscf

 

Download Python script

 

Download Notebook

 

View on GitHub