Source code for pennylane.qchem.matrices

# Copyright 2018-2021 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 the functions needed for computing matrices.
"""
# pylint: disable= too-many-branches
import itertools as it

import autograd.numpy as anp

from .integrals import (
    attraction_integral,
    kinetic_integral,
    moment_integral,
    overlap_integral,
    repulsion_integral,
)


[docs]def mol_density_matrix(n_electron, c): r"""Compute the molecular density matrix. The density matrix :math:`P` is computed from the molecular orbital coefficients :math:`C` as .. math:: P_{\mu \nu} = \sum_{i=1}^{N} C_{\mu i} C_{\nu i}, where :math:`N = N_{electrons} / 2` is the number of occupied orbitals. Note that the total density matrix is the sum of the :math:`\alpha` and :math:`\beta` density matrices, :math:`P = P^{\alpha} + P^{\beta}`. Args: n_electron (integer): number of electrons c (array[array[float]]): molecular orbital coefficients Returns: array[array[float]]: density matrix **Example** >>> c = np.array([[-0.54828771, 1.21848441], [-0.54828771, -1.21848441]]) >>> n_electron = 2 >>> mol_density_matrix(n_electron, c) array([[0.30061941, 0.30061941], [0.30061941, 0.30061941]]) """ p = anp.dot(c[:, : n_electron // 2], anp.conjugate(c[:, : n_electron // 2]).T) return p
[docs]def overlap_matrix(basis_functions): r"""Return a function that computes the overlap matrix for a given set of basis functions. Args: basis_functions (list[~qchem.basis_set.BasisFunction]): basis functions Returns: function: function that computes the overlap matrix **Example** >>> symbols = ['H', 'H'] >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False) >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554], >>> [3.42525091, 0.62391373, 0.1688554]], requires_grad=True) >>> mol = qml.qchem.Molecule(symbols, geometry, alpha=alpha) >>> args = [alpha] >>> overlap_matrix(mol.basis_set)(*args) array([[1.0, 0.7965883009074122], [0.7965883009074122, 1.0]]) """ def overlap(*args): r"""Construct the overlap matrix for a given set of basis functions. Args: args (array[array[float]]): initial values of the differentiable parameters Returns: array[array[float]]: the overlap matrix """ n = len(basis_functions) matrix = anp.eye(n) for (i, a), (j, b) in it.combinations(enumerate(basis_functions), r=2): args_ab = [] if args: args_ab.extend(arg[[i, j]] for arg in args) integral = overlap_integral(a, b, normalize=False)(*args_ab) o = anp.zeros((n, n)) o[i, j] = o[j, i] = 1.0 matrix = matrix + integral * o return matrix return overlap
[docs]def moment_matrix(basis_functions, order, idx): r"""Return a function that computes the multipole moment matrix for a set of basis functions. Args: basis_functions (list[~qchem.basis_set.BasisFunction]): basis functions order (integer): exponent of the position component idx (integer): index determining the dimension of the multipole moment integral Returns: function: function that computes the multipole moment matrix **Example** >>> symbols = ['H', 'H'] >>> geometry = np.array([[0.0, 0.0, 0.0], [2.0, 0.0, 0.0]], requires_grad = False) >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554], >>> [3.42525091, 0.62391373, 0.1688554]], requires_grad=True) >>> mol = qml.qchem.Molecule(symbols, geometry, alpha=alpha) >>> args = [alpha] >>> order, idx = 1, 0 >>> moment_matrix(mol.basis_set, order, idx)(*args) tensor([[0.0, 0.4627777], [0.4627777, 2.0]], requires_grad=True) """ def _moment_matrix(*args): r"""Construct the multipole moment matrix for a given set of basis functions. Args: args (array[array[float]]): initial values of the differentiable parameters Returns: array[array[float]]: the multipole moment matrix """ n = len(basis_functions) matrix = anp.zeros((n, n)) for (i, a), (j, b) in it.combinations_with_replacement(enumerate(basis_functions), r=2): args_ab = [] if args: args_ab.extend(arg[[i, j]] for arg in args) integral = moment_integral(a, b, order, idx, normalize=False)(*args_ab) o = anp.zeros((n, n)) o[i, j] = o[j, i] = 1.0 matrix = matrix + integral * o return matrix return _moment_matrix
[docs]def kinetic_matrix(basis_functions): r"""Return a function that computes the kinetic matrix for a given set of basis functions. Args: basis_functions (list[~qchem.basis_set.BasisFunction]): basis functions Returns: function: function that computes the kinetic matrix **Example** >>> symbols = ['H', 'H'] >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False) >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554], >>> [3.42525091, 0.62391373, 0.1688554]], requires_grad=True) >>> mol = qml.qchem.Molecule(symbols, geometry, alpha=alpha) >>> args = [alpha] >>> kinetic_matrix(mol.basis_set)(*args) array([[0.76003189, 0.38325367], [0.38325367, 0.76003189]]) """ def kinetic(*args): r"""Construct the kinetic matrix for a given set of basis functions. Args: args (array[array[float]]): initial values of the differentiable parameters Returns: array[array[float]]: the kinetic matrix """ n = len(basis_functions) matrix = anp.zeros((n, n)) for (i, a), (j, b) in it.combinations_with_replacement(enumerate(basis_functions), r=2): args_ab = [] if args: args_ab.extend(arg[[i, j]] for arg in args) integral = kinetic_integral(a, b, normalize=False)(*args_ab) o = anp.zeros((n, n)) o[i, j] = o[j, i] = 1.0 matrix = matrix + integral * o return matrix return kinetic
[docs]def attraction_matrix(basis_functions, charges, r): r"""Return a function that computes the electron-nuclear attraction matrix for a given set of basis functions. Args: basis_functions (list[~qchem.basis_set.BasisFunction]): basis functions charges (list[int]): nuclear charges r (array[float]): nuclear positions Returns: function: function that computes the electron-nuclear attraction matrix **Example** >>> symbols = ['H', 'H'] >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False) >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554], >>> [3.42525091, 0.62391373, 0.1688554]], requires_grad=True) >>> mol = qml.qchem.Molecule(symbols, geometry, alpha=alpha) >>> args = [alpha] >>> attraction_matrix(mol.basis_set, mol.nuclear_charges, mol.coordinates)(*args) array([[-2.03852057, -1.60241667], [-1.60241667, -2.03852057]]) """ def attraction(*args): r"""Construct the electron-nuclear attraction matrix for a given set of basis functions. Args: args (array[array[float]]): initial values of the differentiable parameters Returns: array[array[float]]: the electron-nuclear attraction matrix """ n = len(basis_functions) matrix = anp.zeros((n, n)) for (i, a), (j, b) in it.combinations_with_replacement(enumerate(basis_functions), r=2): integral = 0 if args: args_ab = [] if r.requires_grad: args_ab.extend(arg[[i, j]] for arg in args[1:]) else: args_ab.extend(arg[[i, j]] for arg in args) for k, c in enumerate(r): if c.requires_grad: args_ab = [args[0][k]] + args_ab integral = integral - charges[k] * attraction_integral( c, a, b, normalize=False )(*args_ab) if c.requires_grad: args_ab = args_ab[1:] else: for k, c in enumerate(r): integral = ( integral - charges[k] * attraction_integral(c, a, b, normalize=False)() ) o = anp.zeros((n, n)) o[i, j] = o[j, i] = 1.0 matrix = matrix + integral * o return matrix return attraction
[docs]def repulsion_tensor(basis_functions): r"""Return a function that computes the electron repulsion tensor for a given set of basis functions. Args: basis_functions (list[~qchem.basis_set.BasisFunction]): basis functions Returns: function: function that computes the electron repulsion tensor **Example** >>> symbols = ['H', 'H'] >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False) >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554], >>> [3.42525091, 0.62391373, 0.1688554]], requires_grad=True) >>> mol = qml.qchem.Molecule(symbols, geometry, alpha=alpha) >>> args = [alpha] >>> repulsion_tensor(mol.basis_set)(*args) array([[[[0.77460595, 0.56886144], [0.56886144, 0.65017747]], [[0.56886144, 0.45590152], [0.45590152, 0.56886144]]], [[[0.56886144, 0.45590152], [0.45590152, 0.56886144]], [[0.65017747, 0.56886144],[0.56886144, 0.77460595]]]]) """ def repulsion(*args): r"""Construct the electron repulsion tensor for a given set of basis functions. Permutational symmetries are taken from [D.F. Brailsford and G.G. Hall, International Journal of Quantum Chemistry, 1971, 5, 657-668]. Args: args (array[array[float]]): initial values of the differentiable parameters Returns: array[array[float]]: the electron repulsion tensor """ n = len(basis_functions) tensor = anp.zeros((n, n, n, n)) e_calc = anp.full((n, n, n, n), anp.nan) for (i, a), (j, b), (k, c), (l, d) in it.product(enumerate(basis_functions), repeat=4): if anp.isnan(e_calc[(i, j, k, l)]): args_abcd = [] if args: args_abcd.extend(arg[[i, j, k, l]] for arg in args) integral = repulsion_integral(a, b, c, d, normalize=False)(*args_abcd) permutations = [ (i, j, k, l), (k, l, i, j), (j, i, l, k), (l, k, j, i), (j, i, k, l), (l, k, i, j), (i, j, l, k), (k, l, j, i), ] o = anp.zeros((n, n, n, n)) for perm in permutations: o[perm] = 1.0 e_calc[perm] = 1.0 tensor = tensor + integral * o return tensor return repulsion
[docs]def core_matrix(basis_functions, charges, r): r"""Return a function that computes the core matrix for a given set of basis functions. The core matrix is computed as a sum of the kinetic and electron-nuclear attraction matrices. Args: basis_functions (list[~qchem.basis_set.BasisFunction]): basis functions charges (list[int]): nuclear charges r (array[float]): nuclear positions Returns: function: function that computes the core matrix **Example** >>> symbols = ['H', 'H'] >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False) >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554], >>> [3.42525091, 0.62391373, 0.1688554]], requires_grad=True) >>> mol = qml.qchem.Molecule(symbols, geometry, alpha=alpha) >>> args = [alpha] >>> core_matrix(mol.basis_set, mol.nuclear_charges, mol.coordinates)(*args) array([[-1.27848869, -1.21916299], [-1.21916299, -1.27848869]]) """ def core(*args): r"""Construct the core matrix for a given set of basis functions. Args: args (array[array[float]]): initial values of the differentiable parameters Returns: array[array[float]]: the core matrix """ if r.requires_grad: t = kinetic_matrix(basis_functions)(*args[1:]) else: t = kinetic_matrix(basis_functions)(*args) a = attraction_matrix(basis_functions, charges, r)(*args) return t + a return core